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ABSTRACT 

This paper presents accurate absolute positions from a 24 GHz Very Long Baseline Array (VLBA) search 
for compact extragalactic sources in an area where the density of known calibrators with precise coordinates 
is low. The goals were to identify additional sources suitable for use as phase calibrators for galactic sources, 
determine their precise positions, and produce radio images. In order to achieve these goals, we developed a 
new software package, VIAiA , for determining group delays from wide-band data with much lower detection 
limit. With the use of VIAiA we have detected 327 sources out of 487 targets observed in three 24 hour VLBA 
experiments. Among the 327 detected objects, 176 are within 10 degrees of the Galactic plane. This VGaPS 
catalogue of source positions, plots of correlated flux density versus projected baseline length, contour plots, 
as well as weighted CLEAN images and calibrated visibility data in FITS format, are available on the Web 
at |http : / /astrogeo . org/vgaps| Approximately one half of objects from the 24 GHz catalogue were 
observed at dual band 8.6 GHz and 2.3 GHz experiments. Position differences at 24 GHz versus 8.6/2.3 GHz 
for all but two objects on average are strictly within reported uncertainties. We found that for two objects with 
complex structure positions at different frequencies correspond to different components of a source. 
Subject headings: astrometry — catalogues — surveys 



1. INTRODUCTION 

The method of very long baseline interferometry (VLBI), 
first proposed by Matveenko et al. (1965), has numerous ap- 
plications in the areas of high resolution imaging, differen- 
tial astrometry, absolute astrometry, space geodesy, and space 
navigation. Because the turbulence in the neutral atmosphere 
and ionospheric fluctuations set a limit of coherent averaging 
at typically 1 to 10 minutes, depending on frequency, the de- 
tection of weak sources that require longer integration is not 
possible. 

To overcome this limitation, the majority of VLBI obser- 
vations are made in phase referencing mode: the telescopes 
of the array slew rapidly between a weak target and a nearby 
strong calibrator. The phase changes of the calibrator trace 
the fluctuations in the atmosphere, and when they are sub- 
tracted from the phase of the target, the residual phases are 
essentially free from fluctuations caused by the atmosphere, 
and the target integration time can be extended almost indefi- 
nitely, enabling detection and imaging of weak objects. This 
technique is called phase-referencing. 

The technique of phase referencing allows also to deter- 
mine the precise differential position of a target with respect 
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to a calibrator with accuracy reaching 0.05 mas or better, even 
with moderate signal-to-noise ratio (SNR). The advantage of 
differential astrometry over absolute astrometry is that the 
contribution of unaccounted propagation delays and errors in 
station positions is diluted by a factor of the target-to-source 
separation in radians. Either the target or calibrator may be 
observed in a narrow band. 

The ability to image weak sources and determine their po- 
sitions accurately with respect to a nearby calibrat or have 
made phase referencing very popular According to IWrobell 
(|200^. in 2003-2008, 63% of VLBA observations used this 
technique. However, to make phase referencing possible, a 
dense catalogue of phase calibrators is needed such that a suit- 
able calibrator will be found within 2-3° of any target. And 
for precise differential astrometry, the calibrator position ac- 
curacy of a few milliarcsecond is needed. Efforts to create 
such a catalogue of calibrators commenced in t he 1980s under 
the N ASA's Crustal Dynamic Project program ( iRvan & ClarkI 
Il987h th at ultimately res ulted in the ICRF catalogue of 608 
sources dMa et all 119981) . Later, over 6000 sources were 
observed in the framework of the VLBA Calibrator Survey 
(VCS) program ( Beasle y et al.l 120021 iF bmalont et al. 2001 
Petrov et all I2005L 120061 iKova lev et ah 2007, Petro v et alJ 



2007cl) . the VLBA RDV program tPeti'ov et a l. 2009). and the 



continuing Aust ralian Long Baselin e Array Calibrator Survey 
(LCS) program (iPetrov et al.ll201 lb . By the end of 2010, the 
number of known calibrators with position accuracies better 
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than 5 mas surpassed 4600. The probability of finding a cali- 
brator within 2° of any direction is currently 64%, and within 
3° is 90%. However, the distribution of calibrators is not uni- 
form on the sky due to several factors: the large scale structure 
of the Universe; the location of most observing stations in the 
northern hemisphere; and obscuration and confusion within 
5° of the Galactic plane. 

Finding calibrators in the Galactic plane region is espe- 
cially difficult for several reasons. First, the region is filled 
with many galactic objects, and surveys from single antennas 
or km-sized arrays, which are needed to find c alibrator can- 
didate s, avoid this r egion; for example J VAS (IWrobel et alJ 
[T998h and AT20G iMurphv et al.l l20Tol) . Therefore, the 
pool of available candidates is limited near the plane. Sec- 
ond, many potential candidates with flat spectra are extended 
galactic objects, such as planetary nebulae or HII regions, that 
cannot be detected by VLBI. Finally, the apparent angular size 
of extragalactic objects observed through high plasma density 
near the Galactic plane are broadened by Galactic scattering, 
and cannot be detected at low frequencies on baselines longer 
than several thousand kilometers. 

Nevertheless, a dense grid of calibrators in the Galactic 
plane is needed for studies of compact galactic objects in 
both continuum emission (pulsars. X-ray binaries) and line 
emission (egs. water masers, methanol masers. Hydrogen 
absorption). Some extragalactic sources near the Galactic 
plane have b een associated with Ferm/'-detected 7-ray objects 
(lAbdo etal . 2010a b) only through VLBI ca librator surveys as 
it was suggested and successfully shown by lKovalevI (12009^. 

In order to increase the density of calibrators in the Galactic 
plane and to search for suitable calibrators within 2° of known 
water masers, we developed an observing strategy that com- 
bine the VERA and VLBA arrays. First, we systematically 
screened 2462 sources with declinations > -40° using the 4- 
station VLBI array VERA (Honma et al. 2003) at 22 GHz (K- 
band) by observing each object in two 120-second long scans. 
We detected 533 objects , 180 of them new at K-band, and 
these results are given in iPetrov et al.l ( l2007bl) . bmce precise 
determination of parallaxes and proper motions of sources 
with water-maser emission is one of the main targets of the 
VERA project, potential calibrators near known water masers 
were preferentially included in the observations. 

This paper describes the follow-up VLBA observations at 
24 GHz of 487 radio sources in order to determine their pre- 
cise positions and images. We denote these observations as 
the VLBA Galactic Plane Survey (VGaPS), and the results of 
this observing campaign are described in this paper VLBA 
observations, source selection and the scheduling algorithm 
are given in section ^ The data analysis procedure is pre- 
sented in ^ For analysis of these observations we have de- 
veloped a new approach for wide-band fringe sear ch an d as- 
trometric analysis which is described in detail in ^3.21 The 
validation of the results made with the new approach is given 
in ^ The images are described in section ^ the source 
position catalog is listed in ^ and the K- and S/X-band as- 
trometric VLBI positions are compared in ^ The results are 
summarized in section ^ 

2. SOURCE SELECTION AND OBSERVATIONS 

The VLBA observations were made at 24 GHz for several 
reasons. First, Galactic scattering at low frequencies broad- 
ens images and degrades source positions. Even at the fre- 
quency of 8.6 GHz, many calibrators near the Galactic plane 
are too broadened to be useful. Second, many galactic radio 



targets are associated with H2O maser emission at 22.5 GHz. 
Since the correlated flux density of calibrators generally de- 
creases with increasing frequency, a good calibrator at 2.3 and 
8.6 GHz may not be sufficiently compact or bright for use at 
22 GHz. Third, source positions at 22 GHz may not nec- 
essarily be the same as those at 8.6 GHz because of the 
effects of frequency dependent source structure, especially 
for multi-component obj ects, and because of frequency de- 
pendent core shif ts (e.g ., lLobanov|[T998l iKovalev et al1l2008l 
ISokolovskv etaHIIOll . 

We considered the task of extending the pool of calibrators 
for Galactic astrometry more broadly: not only to increase 
the list of compact extragalactic radio sources within 10° of 
the Galactic plane, but also to re-observe known sources at 
K-band that are either in the Galactic plane or close to known 
masers at higher galactic latitudes. Water masers are the main 
target of the VERA project, and determination of their paral- 
laxes and proper motions is important for studying the three- 
dimensional structure and dynamics of the Galaxy's disk and 
bulge, and for revealing the true shape of the bulge and spi- 
ral arms, its precise rotation curve and the distribution of dark 
matter. Dual-beam K-band VERA observations require cal- 
ibrators with positions known to the milliarcsec level within 
2.° 2 of the target. 

2.1. Source selection 

VLBI can detect emission only from a compact component 
of a source, which should be bright enough to be detected 
above the noise background. Chances to find a radio source 
sufficiently compact for VLBI detection are significantly en- 
hanced if information about source spectra is available. For a 
majority of sources with flat spectrum, i.e. with the spectral 
index a > -0.5 (S ~ /"), synchrotron emission from a com- 
pact core often dominates, so such sources are often compact. 
In regions more than 10° from the Galactic plane, compari- 
son of surveys at different frequencies can be used to select 
the sources with flat spectra. However, in the crowded Galac- 
tic plane, source misidentification and confusion from surveys 
at different frequencies often result in ambiguous informa- 
tion about a source spectrum. Also, there are many extended 
sources with thermal emission, such as planetary nebulae and 
HII regions, that have flat spectrum but are too extended for 
VLBI detection. 

To determine if a candidate source was sufficiently strong 
and compact to b e a VLBI calibrator, we first used the VERA 
array at 22 GHz ( IPetrov et al.ll2007bl) . In the survey --1400 
objects within 6° of the Galactic plane and ~'1000 objects 
within 2° of known maser sources outside the Galactic plane 
were observed with the four-element VERA array at baselines 
of 1000-2000 km long in two scans each, and approximately 
20% were detected. Among the 533 detected sources, 305 
objects are with galactic latitude \b\ < 10°, and 228 objects 
are with \b\ > 10°. The list of detected sources included all 
objects with the probability of false detection < 0.1. These 
objects formed the first set of candidate calibrators. 

We also added 239 objects in the Galactic plane and 44 
objects outside Galactic plane, selected on the basis of their 
spectra using the Astrophysical CATalogs support System 
CATS ( Verkhod anov et al. 2005). This database currently in- 
cludes catalogues from 395 radio astronomy surveys. We se- 
lected all entries with sources within a 20"-radius and with 
measurements of flux density at least at three frequencies in 
the range of 1.4-100 GHz. We fit a straight line to the log- 
arithm of the spectrum and then estimated the spectral index 
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and the flux density extrapolated to 24 GHz. We selected 451 
objects with extrapolated correlated flux densities at K-band 
> 200 mJy, spectral index flatter than -0.5 and \b\ < 10°. 
Then, we scrutinized this list and removed sources within 30" 
of known planetary nebulae or HII regions and sources with 
anomalous spectra that indicated a possible misidentification. 
The remaining list was added to the pool of calibrators. This 
list contained also Sagittarius A* since it has never before 
been observed with VLBI in the absolute astrometry mode. 

This set of sources formed the pool of 8 16 candidate objects 
to be followed-up with the VLBA. Because of the large num- 
ber of candidate objects, a priority from 1 to 4 was assigned to 
each source. The first priority was given to 180 new sources 
with correlated flux densities in the range 100-300 mJy de- 
tected with VERA that have never been observed with VLBI 
before. The second was given to sources in the Galactic plane 
detected with VERA; the third was given to other sources de- 
tected with VERA; and the least priority was given to sources 
outside the Galactic plane not detected with VERA. All new 
sources detected with VERA were scheduled. Source out- 
side of the Galactic plane were included for logistical reasons. 
Since the Galactic plane has an inclination with respect to the 
equator of 62°.6, the density of candidate sources near the 
Galactic plane is a non-uniform function of right ascensions. 
We thus included sources outside the Galactic plane in the 
right ascension regions that had significantly fewer candidate 
sources than others in order to avoid gaps in the schedules. 

In addition to the ta rget sources, we a lso selected 56 objects 
from the K/Q survey (iLanyi et al.ll2010h to serve as amplitude 
and atmosphere calibrators. These are the brightest sources at 
K-band with precisely known positions and publicly available 
images in FITS-format. 

2.2. Scheduling algorithm 

A sufficiently good preliminary VLBA observing schedule 
for three 24-hour sessions was prepared automatically. Then 
it was manually adjusted in order to produce a more efficient 
final schedule. 

The three sets of information needed in order to compile 
the schedule were: 1) the list of the 816 target sources with 
their J2000 positions and priority level; 2) the list of the 56 
calibration sources at 24 GHz; and 3) the two or three side- 
real times ranges (at the array center at station PIETOWN) that 
each source should be observed. These times were a function 
of the source declination. For example, a source with S > 50° 
could be observed three times with very flexible time ranges; 
whereas, a source with S < -35°, must have tight ranges in or- 
der to be observed by at least 8 VLBA antennas at an elevation 
higher than ten degrees. 

The scheduling goal was to observe each source for three 
scans of 120 seconds length, unless it was south of declina- 
tion -25°, in which case only two scans were scheduled. With 
an average overhead of about 45 sec between sources, about 
1550 scans over the three days (72 hours) could be sched- 
uled. Every 90 minutes, four scans were reserved for atmo- 
sphere calibrators. The algorithm then began filling in observ- 
ing slots, taking the highest priority sources first, and those at 
the lowest declinations, since these have minimal scheduling 
flexibility. The algorithm scheduled sources that were rela- 
tively close in the sky in order to minimize slew times, which 
can be as long as three minutes. 

The slots containing the calibrators were chosen in order 
to maximize the range of elevations for the VLBA antennas. 
In practice, this meant maximizing elevations of observed 



Table 1 

VLBA Target Sources from VERA Observations 



group 


pool 


scheduled 


detected 


ratio 


Galactic, VERA 


305 


184 


140 


76% 


Non-galactic, VERA 


228 


151 


133 


88% 


Galactic others 


239 


108 


36 


33% 


Non-galactic, others 


44 


44 


15 


34% 


Total 


816 


487 


327 


67% 



sources at SC-VLBA, MK-VLBA and PIETOWN. It was impor- 
tant to have one low elevation observation for all telescopes, 
and this could be scheduled by observing a source either in the 
far north and/or in the far south. All three days were sched- 
uled at the same time, since it did not matter on which day 
any of the three sidereal time slots occurred for a source. 

The final schedule optimization was done 'by hand' and 
consisted of several steps. First, some sources could be placed 
in only one or two slots and would have to be removed if an- 
other slot could not be found. Since about 10% of the slots 
could not be filled using the automatic algorithm because the 
source list was not uniformly distributed over the sky, addi- 
tional slots were usually found to complete a source's sched- 
ule requirement. This often meant bending some of the rules 
or moving calibration sources or blocks by five to fifteen min- 
utes. In regions where the target source density was large, 
scan integrations were decreased from 120 sec to 1 10 sec. All 
sources with priority 1 were included. 

The next stage insured that the calibration scans were opti- 
mized in order to obtain good elevation coverage for the an- 
tennas. The purpose of these observations was 1) to serve as 
amplitude and bandpass calibrators; 2) to improve robustness 
of estimates of the path delay in the neutral atmosphere; and 
3) to tie the source positions of new s ources to existin g cat- 
alogues, such as the ICRF catalogue (iMa et al.iri998h . The 
final stage of optimization tweaked the observing schedule in 
obvious ways in order to save slewing time since the schedul- 
ing algorithm did not minimize slew times, and because of 
the above adjustments to insure proper observing coverage for 
each source. The schedule for each session was then carefully 
checked using the NRAO SCHED program to insure that each 
scan had sufficient integration time, and no more than one 
of the antennas was below ten degree elevation (except for 
sources south of -35°). 

The results of the scheduling and ultimate detection for 
each priority group is given in Table [T] The number of tar- 
get sources selected in each group is given in column 2, those 
scheduled in column 3, those detected by the VLBA in col- 
umn 4, and the detection rate in the last column. The first two 
groups are sources detected in VERA Fringe Search Survey, 
the last two groups contain other objects. The table does not 
include the 56 sources used as atmosphere calibrators. 

2.3. Observations 

The VGaPS observations were carried out in three 24 hour 
observing sessions at the VLBA on 2006 June 04, 2006 June 
11, and 2006 October 20. Each target source was observed 
in several scans: 3 sources in 1 scan, 356 sources in 2 scans, 
1 source in 3 scans, 124 in 4 scans and 3 sources in 5 scans. 
The scan durations were 100-120 seconds. In total, anten- 
nas spent 57% of time on target sources. In addition to these 
target objects, 56 strong sources previously observed at K- 
band were taken from the astrometric and geodetic catalogue 
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Table 2 

The lower edge IF frequencies (in GHz) 



IFl 24.20957 

IF2 24.22257 

IF3 24.26157 

IF4 24.33957 

IF5 24.48257 

IF6 24.58657 

IF7 24.65157 

IF8 24.67757 



2004f_astrcC]. 

The data digitized at four levels were recorded with the rate 
of 256 Mbit/s in eight 8 MHz wide intermediate frequencies 
(IF) bands spread over a bandwidth of 476 MHz (Table |2]i. 
The frequencies were selected to minimize the amplitude of 
sidelobes of the Fourier transform of the bandpass. 

3. DATA ANALYSIS 

The data were correlated in Socorro on the VLBA hardware 
correlator The correlator output contains the complex spectra 
of the autocorrelation function and the spectrum of the cross- 
correlation function for each accumulation period. The accu- 
mulation period was chosen to be 0. 1 3 1 072 s, and the spectral 
resolution was set to 0.5 MHz, i.e. 16 spectral channels per 
IF 

After correlation, the data were stored in a file compliant 
with the FITS-IDI specifications (Eric Greisen, NRAO memo 
Nll43). Small a priori amplitude corrections were applied; 
bad data were flagged (sometimes bad data are not found until 
later processing); and the definitions of the reference frequen- 
cies were modified. 

Further data analysis involves computation of group delays 
for each scan and each baseline, computation of theoretical 
path delay, and then fitting parameters of the linear model into 
the differences between the observed and theoretical delays 
using least squares (LSQ). 

3.1. Traditional narrow-band fringe fitting algorithm 

The data set then contains 128 data streams for each of the 
45 baselines if all ten VLBA antennas are operating: 8 IF's, 
each with 16 equally-space frequency channels. After the 
complex bandpass calibration, the phase difference among all 
of the frequency streams remains unchanged, and the relevant 
residual phase parameters associated with any scan are 1) the 
residual phase at the scan midpoint, 2) the average group de- 
lay (phase gradient with frequency), and 3) the average delay 
rate (delay gradient with time), for each antenna. These pa- 
rameters are called the residual phase terms. These parame- 
ters are estimated with a fringe fitting procedure. Generally, 
one antenna is chosen as the reference — a well-behaved one 
near the array center — so that set of a residual phase terms 
are associated for each scan and all other antennas. These pa- 
rameters are functions of many astrometric quantities (source 
position, site positions, antenna-based tropospheric path de- 
lays. Earth Orientation parameters, etc), which can be deter- 
mined from analysis systems like Calc/Solve from data ob- 
tained from carefully prepared observing schedules. Source 
structure and dispersive phase effects (e.g. ionosphere) pro- 
duce a non-linear phase/frequency relationship. 



The algorithm that determines the residual phase terms 
is implemented in the AIPS task FRJNG, and in the past, 
all VLBA observations under the absolute astrometry and 
geodesy programs, such as VCS, RDV, and K/Q , were pro- 
cessed with use of this software jGreisen 1 [19881) . WhenlFs 
are spread over a wide band with gaps, the AIPS algorithm de- 
termines the residual phase in two steps: first, for each of the 8 
IF's the residual phase, single-band group delay and rates are 
independently obtained. Then, the residual phases of each IF 
are fit to a linear phase versus frequency term to produce the 
group delay, using A IPS program MBDL Y. This procedure is 
described in detail in lPetrov et alj ( 120091) . 



3.2. Wide-band baseline-based fringe fitting algorithm 

The two-step approach has a substantial shortcoming: it re- 
quires fringe detection for each individual IF independently 
within a narrow band. Using all IFs simultaneously for co- 
herent averaging, we can detected a source with the amplitude 
VA^ smaller than using only one IF. This degradation of the 
detection limit does not pose a problem for most geodesy ob- 
servations, or for absolute astrometry of bright sources, since 
the target SNR is usually very high, but it impacts signifi- 
cantly absolute astrometric experiments of weak sources. The 
traditional AIPS algorithm did not detect a sizeable fraction 
of the sources observed in the VGaPS experiment. 

This limitation motivated us to develop an advanced algo- 
rithm for wide-band fringe fitting across all of the IF's within 
the band. For logistical reasons, instead of augmenting AIPS 
with the new task, we decided to develop from scratches 
the new software package called T^MA □ that is supposed 
to replace the AIPS for processing absolute astronomy and 
geodesy experiments. Here we outline the method of wide- 
band fringe search used for processing this experiment. 



3.2.1. Spectrum re-normalization 

Digitization of the input signal and its processing with 
a digital correlator causes an amplitude distortion with re- 
spect to an ideal analogue system. As documented by Koga^ 
many effects distort both cross-correlation and auto- 
correlation spectrum exactly the same way. Therefore, if we 
divide the cross-correlation spectrum by the auto-correlation 
spectrum averaged over time and over frequencies within each 
IF, we will remove these distortions. However, there are two 
effects that affect cross- and auto-spectrum differently be- 
cause the amplitude of the cross-spectrum is very low, and 
the amplitude of the auto-correlation spectrum is close to 1 . 

The non-linear amplitu de distortion o f the digitized signal 
was studied in depth by iKoganI (Il998h who derived a gen- 
eral expression for the correlation coefficient of the digitized 
signal as a function of the correlation coefficient of the hypo- 
thetical analogue signal. In the absence of fringe stopping, the 
output correlation coefficient po,,, is expressed via the correla- 



*'http : / / ast rogeo ■ org/ rf c| 

2 Available at ftp : / / ftp ■ aoc ■ nrao ■ edu/pub/ software / Xaips/SSMEKlBlgi^Byt^l-ygMEftHdlgte.Bfib/pimal 
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tion coefficient for an analogue case p as 



Pout — 2 k 



1 



■.dp 



+ 2k(«-1) 



+ K{n-\f 



_pj-2pp,i.2+v5 _V|+2pv,r2+ _ 



where k is the normalization coefficient, « is the ratio of two 
level of quantization, vj and \'2 are the dimensionless digitizer 
levels in units of variance of the inpu t signal. Their numeri- 
cal values compiled from the paper of lKoganI (119931) are pre- 
sented in table [3] 



Table 3 

Numerical coefficients in integral[T]for three cases of the number of bits per 
sample: (1,1), (1,2), (2,2) 





n 




V2 


K 


(1,1) 


1.0 


0.0 


0.0 


0.3803 


(1,2) 


3.3359 


0.0 


0.9816 


0.05415 


(2,2) 


3.3359 


0.9816 


0.9816 


0.07394 



When p ^ 1, the dependence Pomip) becomes linear: pom = 
2/tt ■ p 0.6366 p for the case of 1-bit sampling, and pom ~ 
0.8825 p for the case of 2-bit sampling. Therefore, distortion 
of the cross-spectrum is proportional to Pout/p, and dividing 
the cross-spectrum by 0.6366 or 0.8825 we eliminate the dis- 
tortion of the fringe amplitude introduced by the digitization. 

However, the amplitude of the autocorrelation spectrum is 
close to 1, and we cannot ignore non-linearity of PouAp)- To 
correct the autocorrelation spectrum fo r digitization distor- 
tion, we follow the procedure outlined by lKogaril(ll995h . First, 
the auto-correlation is inverse Fourier transformed. It should 
be noted that the correlator provides the autocorrelation for 

spectral channels for non-negative frequencies from to 
N-l. We restore the autocorrelation at the A^-th channel by 
linear extrapolation using the A^- 2 -th and A^- 1 -th values of 
the spectrum and set the spectrum for A^- 1 negative frequen- 
cies to zero. The dimension of the Fourier transform is 2A^. 
Second, the result of the transformation, the auto-correlation 
coefficient versus time lag, is de-tapered, i.e. divided by the 
self-convolution of the weighting function. The VLBA corre- 
lator normally uses uniform weighting, i.e. uses weight 1 for 
all points. The self-convolution of the uniform weighting is a 
triangle function A(i): 



A(/): 




if |/|<A^ 
otherwise 



Third, the correlation function is divided by its maximum 
that is found at the zero lag. Fourth, each correlation coef- 
ficient is divided by Pomip)- Fifth, the correlation function 
is multiplied back by the stored value at the zero lag. Sixth, 
the correlation function is again tapered by multiplying it by 
A(i). Finally, we perform the Fourier transform of the cor- 
rected correlation function and get the re-normalized auto- 
spectrum, free from digitization distortion. The square root 
of the product of the auto-correlation spectra from two sta- 
tions of a baseline, averaged over time and frequency within 



each IF, gives us an estimate of the fringe amplitude scale fac- 
tor But before dividing the cross-correlation spectrum by this 
scale-factor we have to take into account a specific effect of 
the hardware VLBA correlatofl An insufficient number of 
bits in internal correlator registers resulted in a decrease of 
the amplitude when its is large, and when i t happens, the au- 
tocorrelation spectra are corrupted. iKoganI ( 11993 ') suggested 
the following model for accounting for this effect: 



F = l + - 



w 



ASAV, 



(2) 



where w is the weight of the spectrum data equal to the ratio 
of processed samples to the total number of samples in an ac- 
cumulation period, A is the accumulation period length, S is 
2 when the correlator processed single polarization data and 
1 if both polarizations were correlated, and Vs is the visibility 
scale factor provided by the correlator. We divide the autocor- 
relation spectrum by the F factor. 

3.2.2. Coarse fringe search 

The correlator output provides autocorrelation and cross- 
correlation spectra for each pair of baselines and each scan. 
The cross-correlation spectrum is computed at a uniform two- 
dimensional grid of accumulation periods and frequencies and 
is accompanied with weights that are the ratio of the number 
of processed samples in each accumulation period to the nom- 
inal number of samples. 

The fringe fitting procedure searches for phase delay Tp, 
phase delay rate Tp, group delay Tg, and its time derivative fg 
that correct their a priori values used by the correlator model 
in such a way that the coherent sum of weighted complex 
cross-correlation samples Ckj 

k j (3) 

^liuJaTp + i^o-rp((i:-/o) + (u]j-uJo)'rg + (uij-uJo)Tg(h-to)) 

reaches the maximum amplitude. Index k runs over time, and 
index j runs over frequencies, and fo denote angular refer- 
ence frequency within the band and the reference time within 
a scan and Wkj are weights. Function C{Tp,Tg,Tp,T'g) is essen- 
tially non-linear, and we need a really good starting value in 
order to find the global maximum by traditional optimization 
algorithms. We can notice that term InojoTp in expression[3] 
does not depend on the summation indexes, and Tg is usually 
small. Therefore, for the purpose of coarse fringe search we 
simplify expression[3]to 



CiTp,Tg,Tp)e 



-llTTUJoT,, 



k j 



JiuioTpitk-'o) + (Wj-wo)-rj) 



(4) 



For the search of the maximum, the trial functions 
C{Tp,Tg,Tp) are computed on a dense grid of the search space 
Tg,Tp. It follows immediately from expression |4] that \C\ = 
\J-{ctj)\ , where T denotes the two dimensional Fourier trans- 
form. 

The first step of the coarse fringe search is to compute the 
two-dimensional Fast Fourier Transform (FFT) of the matrix 
of the cross-correlation spectrum. The first dimension of the 
matrix runs over time, the second dimension runs over fre- 
quency. The sampling intervals are At//3 and A//7 where 

^ The new VLBA software DiFX correlator does not have this problem. 
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At and A/ are the duration of the accumulation period and 
the spectral resolution respectively, and /? and 7 are integer 
oversampling factors. The elements of the matrix which do 
not have measurements or have discarded measurements are 
padded with zeroes. The dimensions of the matrix are chosen 
to have the power of 2 for gaining the maximum performance 
of the FFT numerical algorithm. 

The oversampling factors greater than 1 are used for miti- 
gation of amplitude losses. The FFT produces the estimates 
of ICI at a discrete grid of group delays and delay rates. If 
the maximum of the amplitude of the coherent sum of the 
cross correlation function samples falls at group delays and 
delay rates between the nodes of the grid, its magnitude will 
be greater than the amplitude of ICI at the nearest grid point 
by a factor L. 

The amplitude loss factor L at the coarse search matrix is 
the integral average over time and frequency: 



1 r'' f k\ 

L=— I cos27r woTp-— — fc/fx 

J-/,/2 V Pt, J 

1 ff"'^ ( I \ 

■J- / cos27r uJoTg-— fdf, 

Jb J-f,l2 \ lib J 



(5) 



where k and / are indexes of the nearest grid nodes, ts is 
the scan duration, is the total bandwidth. The integral 
|5] is easily evaluated analytically, and the maximum losses 
L = sine (7r/(2/?)) • sine (tt/ (27)) are achieved when r'p and 
Tg happen to be just between grid nodes. In the case when 
the oversampling factor is 1, L = 4/7r^ « 0.405. That loss fac- 
tor effectively raises the detection limit by 1/L = 2.467 in the 
worst case. We used the grid 4096 x 4096 which corresponds 
to /3 = 4, 7 = 4 for = 120 s, f,, = 476 MHz. Therefore, the 
maximum loss factor for our experiment is 0.959, which re- 
sults in raising the detection limit by no more than 4.1%. 

The group delay and delay rate that correspond to the maxi- 
mum of the discrete Fourier transform of the cross-correlation 
matrix provide the coarse estimates of group delay and delay 
rate. Their accuracy depends on the grid resolution. The next 
step is to refine their estimates. The first stage of the fine 
search is an iterative procedure that computes ICI at a pro- 
gressively finer three-dimensional grid at the close vicinity of 
the maximum using the discrete 3D Fourier transform. The 
third dimension is group delay rate, omitted during the coarse 
search. Dimensions of the transform for group delay, phase 
delay rate, group delay rate are 3, 3, and 9 respectively. At the 
first step of iterations, the grid runs over ±1 element of the 
coarse grid for group delays and phase delay rates and in the 
range ±2 • 10~" for group delay rate. After each step of itera- 
tions, the grid centered around the maximum element shrinks 
its step by 2. In total, 8 iterations are run. The phase of the 
coherent sum of cross correlation function determined with 
the iterative procedure according to expression|4]is the fringe 
phase with the opposite sign. 

3.2.3. The probability of false detection 

The significance of the fringe amplitude depends on the 
level of noise. In the absence of signal, the amplitude of the 
cross correlation function a has the Rayleigh distribution: 



p{a) = ■ 



(6) 



where a is the standard deviation of the real and imaginary 
part of the cross correlation function and n is the total number 



of spectrum points. In the case when all points of the spec- 
trum are statistically independent, the cumul ative distribution 
function of the coherent sum over « points is ([Thompson et al.l 
I2OOI 



P(a)= l-e V2„ 



(7) 



Differentiating this expression, we find the probability den- 
sity function of the ratio of the amplitude of the coherent sum 
of spectrum samples of the noise to its standard deviation as 



(8) 



p{a) = n^e 2^ \ l-e 2^ 



Under an assumption that all samples are statistically inde- 
pendent, the variance of the noise of the coherent sum of 
samples is scaled as cr = ctj/ -y/fj S,-, where Sr is the sampling 
rate of the recorded signal (6.4 • 10^ in our experiment), is 
the scan duration, and is tr, the variance of an individual sam- 
ple , 1 for a perfect system. 

However, the assumption of statistical independence is an 
idealization. The presence of systematic phase errors, devi- 
ation of the bandpass from the rectangular shape and other 
factors distort the distribution. The deviation from the statis- 
tical independence is difficult to assess analytically. 

We evaluate the variance of the noise by estimating the vari- 
ance over a sample of 32768 random points at the region of the 
Fourier transform of the coherent sum of the cross-correlation 
that does not contain the signal. The indexes of grid elements 
of the sample are produced by using the random number gen- 
erator The sample of amplitudes is ordered, and one half of 
the points greater than the median is rejected. The variance 
over remaining points is computed and an iterative procedure 
is launched that adds back previously rejected points in the 
ascending order of their amplitudes and it updates the mean 
value and variance. The iterations are run till the maximum 
amplitude of the next sample reaches 3.5crfl. The initial re- 
jection and consecutive restoration of points with amplitudes 
> 3.5(Ta ensures that no points with the signal from the source 
affect computation of <Ta and<fl>. The rejection of the tail 
of the amplitude distribution causes a bias in estimate of the 
mean, but the magnitude of the bias is only -2 • 10"'*, which is 
negligible. We define a signal to noise ratio as SNR = a/ <a>. 
It follows immediately from expression|6]that <a>= 

We assume that the a posteriori distribution of the signal to 
noise ratios s = a/aa can be approximated as a function like 
expression[8]with effective parameters fXeff and «eff: 



2 «eff 

p(s) = se 

TT CTeff 



(,-.-4) 



Icff-l 



(9) 



These parameters Ceff and «eff can be found by fitting the 
left tail of the empirical distribution of the signal to noise ra- 
tios. Using their estimates, we can evaluate the probability of 
finding the amplitude less than a if no signal is present, i.e. 
the probability of false detection: 



Pf(s)=l- I p(s)ds=l-— ( 1- 

feff 



(10) 



The low end portion of the empirical SNR distribution and 
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2.0 




SNR 

Figure 1. The low end of the empirical distribution of the achieved SNR for 
the fringe amplitude from results of fringe fitting VLBA data (tilled circles) 
and the fitted curve (thin line) of the theoretical distribution for the case of no 
signal. 

the table of the probability of false detection for VGaPS ex- 
periments are shown in Figure [T] and Table |4] 

Table 4 

The probability of false detection as a function of SNR 



SNR 


Pf{s) 


4.96 


0.3 


5.19 


0.1 


5.61 


0.01 


5.99 


0.001 


6.68 


10-5 



3.2.4. Fine fringe search 

Finally, the group delay, phase delay rate, group delay rate 
and fringe phase at the reference frequency are estimated us- 
ing LSQ in the vicinity of the maximum provided its ampli- 
tude exceeds the detection limit. The goal of the LSQ re- 
finement is to get realistic estimates of statistical errors of fit- 
ted parameters and to account for possible systematic errors 
by analyzing residuals. All cross-correlation spectrum data 
points of a given observation are used in a single LSQ solu- 
tion with weights reciprocal to their variance. 

Determining the variance of the fringe phase of an individ- 
ual point is trivial only in two extreme cases: when SNR <C 1, 
and when SNR ^ L In the first case the fringe phase distri- 
bution becomes uniform in the range of [0,27r] with the vari- 
ance TT / V3 w L 8 1 3 . In the second case expanding the expres- 
sion for fringe phase in the presence of noise </> = arctan 
where Si and S,- are the real and imaginary parts of the signal 
and «, and «, are the real and imaginary parts of the noise, 
into the Taylor series, neglecting terms 0{n/S^), and evaluat- 
ing the variance of the expansion, we get = y^sNR ■ ^'^^ 
the general case, the problem becomes more difficult, since 
the (70 depends on the variance of the noise and on the ampli- 
tude of the signal non-linearly. An analytical solution requires 
evaluation of complicated integrals that are not expressed via 
elementary functions. 

We used the Monte Carlo approach for computing these 
variances. Let us consider a random complex set s = A + 
rir + irii, where A is the amplitude of the simulated signal and 
are independent random variables with Gaussian dis- 




SNR=5 



°° 2 4 6 8 10 

Normalized amplitude A/a 

Figure 2. The variance of fringe phase as a function of the normalized am- 
plitude A I (T„ in the presence of a signal with given SNR. 

tribution. Their variance a was selected in such a way that 
A = ^/^(J SNR. Then for a given SNR, we can compute the 
variance of phase of s as a function of the normalized ampli- 
tude \s\/(j. This is done by generation a long series (1 billion 
points) of the simulated complex signal for a given SNR, com- 
puting the amplitude and phase of the time series, splitting the 
signal into a uniform grid of 128 bins over normalized ampli- 
tude that spans the interval [SNR - 4.5, SNR + 4.5], comput- 
ing the variance of the phase of the simulated signal over all 
points that fall into each bin, and approximating the depen- 
dence of a^{\s\/(T) with a smoothing B-spline of the 3rd order 
over 6 nodes. We computed cr^dil/cr) for SNRs in the range 
[0, 12.7] with steps of 0.1. Examples of this dependence for 
several SNRs are shown in figure |2l The set of B-spline co- 
efficients forms a two dimensional table with axes SNR and 
normalized amplitude that allows evaluating cr^ for a given 
SNR and a given fringe amplitude. 

It should be noted that the SNR of the coarse search is 
related to the amplitude coherently averaged over all valid 
cross-correlation spectrum samples, Wkj, where A', is 

the number of accumulation periods and A^, is the total num- 
ber of spectral channels. The SNR of an individual accumu- 
lation period (the elementary SNR) is ^ J^'j' ^kj times 
less. For our experiments, the typical reduction of the SNR 
is a factor of 340. This means that for almost all sources the 
elementary SNR will be less than 1 . As we have seen previ- 
ously, the distribution of fringe phases at very low SNR's is 
close to uniform with a variance of tt/ \/3. The phase becomes 
uncertain due to the 2tt ambiguity, and the LSQ estimation 
technique loses its diagnostic power. Therefore, the cross- 
correlation function with applied phase, phase delay rate and 
group delays computed during the coarse fringe search have to 
be coherently averaged over time and frequency in segments 
large enough to have sufficiently high SNR over a segment to 
provide an ambiguous phase. As it is seen in Figure [3] the 
fringe phase distribution at SNR=1 is already sharp enough 
for that. Therefore, the number of spectral channels and the 
number of accumulation periods within a segment is chosen 
in such a way that the SNR be at least 1 . Marginally detected 
scans with SNR=5 have 24 segments that average all spectral 
channels within an IF and over 1/3 of the scan interval. 

Using all segments, we determine four fitting parameters p 
using the LSQ; 

p = (A^WAr^A^W<ps, (11) 
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Figure 3. The probability density distribution of fringe phase with SNR=1 
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where A is the matrix of observations, W is the diagonal 
weight matrix and 4>s is the vector phases of the cross- 
correlation function averaged over segments. 

The mathematical expectation of the square of the weighted 
sum of residuals R in the presence of noise e is 

E(R) = Tr(yVCove) - Tr (CoveWA(A^W Ar'A^W) , (12) 

which is reduced to n - m where « is the number of equations 
and m is the number of estimated parameters if the weight 
matrix W is chosen to be (Cove)"' . 

The presence of additive errors, for example fluctuations in 
the atmosphere, will increase Cove and our estimate of the 
error variance based on the amplitude of the spectrum sam- 
ple without knowledge of the scatter of the cross-spectrum 
phases is incomplete. A small additive noise with the vari- 
ance k times less than the amplitude of the signal affects the 
amplitude as 0{k^), but it affects phase as 0{k). One of the 
measures of the incompleteness of the error model is the ratio 
of the square of the weighted sum of residuals to its mathe- 
matical expectation. 

We can extend our error model assuming that the weight 
matrix is >V = Cove"' -q'l, where / is the unit matrix of the 
same dimensions as W and q is the parameter This model is 
equivalent to an assumption that the used squares of weight of 
each segment were less than the true one by some parameter 
q equal to all segments. The mathematical expectation of R 
for this additive weighting model is 

EiR) = n-m+q^[TT(W)-TT{W{AW^A)-^A^)] . (13) 

Inverting [T3l we find q for a given £(/?): 



E{R)-{n-m) 



Tr(W) - Tr(>V(AW2A)"'A 



(14) 



Replacing the mathematical expectation of R with its value 
evaluated from the residuals, we can find the re-weighting pa- 
rameter q for a given solution. Several iterations provides 



Applying an additive 



^ E(R) 

a quick convergence of — 5— to 1. 

R 

reweighting constant results in an increase in estimates of pa- 
rameter uncertainties. It may happen that q becomes imagi- 
nary. This means that W was overestimated. In our analysis 
we set ^ = in that happens. 

3.2.5. Complex bandpass calibration 



Coherent averaging of the cross-correlation spectrum over 
frequencies assumes that the data acquisition system does not 
introduce a distortion of the recorded signal, but this is usu- 
ally not the case. Each intermediate frequency has its own 
arbitrary phase offset and group delay that may vary with 
time. The imperfection of baseband filters results in a non- 
rectangular shape of the amplitude response. 

For calibrating these effects, a rail of narrow-band phase 
calibrator signals with spacing of 1 MHz was injected near the 
receivers. Two tones per IF are extracted by the data acqui- 
sition hardware, and their phase and amplitude are available 
for data analysis. When the phase of the phase calibration sig- 
nal is subtracted from fringe phases, the result is referred to 
the point of injection of the phase calibration signal and this 
procedure is supposed to take into account any phase changes 
that occurred in the signal passing through the data acquisi- 
tion terminal. However, we should be aware of several com- 
plications that emerge when we try to use the benefits of the 
phase calibration signal. First, the phases of two tones of the 
phase calibration signal separated by 6 MHz, as in our ses- 
sions, have ambiguities. Since the instrumental group delay 
may reach several phase turns over a 8 MHz IF, the second 
phase calibration tone is useless without resolving the ambi- 
guity. Second, the phase calibration itself may have a phase 
offset or may become unstable if its amplitude is not carefully 
tuned. Therefore, we need to re-calibrate the phase calibration 
signal itself in order to successfully apply it to data. 

We compute the complex bandpass function for each sta- 
tion, except the reference station, that describes the residual 
instrumental complex bandpass after applying the phase cal- 
ibration signal. The cross-correlation spectrum needs to be 
divided by the complex bandpass in order to correct the in- 
strumental frequency-dependent delay and fading of the am- 
plitude. The procedure for evaluating the complex bandpass 
has several steps. 

First, all data are processed applying the first tone of the 
phase calibration, i.e. the phase of the phase calibration signal 
is subtracted from each phase of the cross-correlation signal 
at a given IF. 

Second, the bandpass reference station is chosen. Then for 
each station, we find an observation at a baseline with the 
reference station that provided the maximum SNR during the 
first run. Then for each IF we average the residual spectrum 
over time and perform a linear fit to the residual phases for 
determination of the instrumental group delay. Using this in- 
strumental group delay, we extrapolate the phase of the phase 
calibration signal of the first phase calibration tone to the fre- 
quency of the second phase calibration tone and resolve its 
phase ambiguity. After that, we re-run the fringe search pro- 
cedure for these scans by applying the phase calibration phase 
to each IF in the form of a linear function of phase versus 
frequency that is computed from two phase calibration tones 
with resolved phase ambiguity. The result of the new fringe 
search gives a new residual spectrum. Then we averaged the 
residual spectrum, over time and overM segments within each 
IF (M = 2 in our experiment). The amplitude of the spectrum 
is normalized by dividing the average amplitude over all IFs. 
The phase and the amplitude of the residual averaged spec- 
trum are approximated with a Legendre polynomial of degree 
5. The result of this approximation as a complex function 
of frequency defines the so-called initial complex bandpass. 
Analysis of residuals of rejected observations helps to diag- 
nose malfunction of the equipment. For example, one or more 
video-converters may fail, which may result in a loss of co- 
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herence. In that case, the part of the affected cross-correlation 
spectrum is masked out. 

Third, we refine the complex bandpass. We select more 
observations with the highest SNR from the first run at all 
baselines (A^=16 in our experiment) and re-run the fringe fit- 
ting procedure with applied phase (but not amplitude!) of the 
initial bandpass. We compute a residual spectrum averaged 
over time and M segments for each processed observation and 
normalize its amplitude. Then we fit a set of 6 coefficients of 
the Legendre polynomial for each station, except the refer- 
ence one, for each IF, for both amplitude and phase to the 
phase and amplitude of the residual spectrum using a single 
LSQ solution. Then the residuals are computed and the obser- 
vations with the maximum absolute value of residual phases 
and residual amplitudes are found separately. If the maxi- 
mum absolute value of residual phase or residual amplitude 
exceeds the predefined limit, the affecting observation is re- 
moved and the solution is repeated. Iterations are run until 
either the absolute values of all remaining residuals are less 
than the predefined limit or the number of observations at a 
given baseline drops below N /2. The fitted Legendre poly- 
nomial coefficients are added to the coefficients of the ini- 
tial bandpass and the result defines the so-called fine complex 
bandpass. 

Fourth, all observations are reprocessed with the fine band- 
pass applied: the phase of the fine complex bandpass of the 
remote station of a baseline is subtracted and the phase of the 
bandpass of the reference station of a baseline is added be- 
fore fringe fitting, and the amplitude is divided by the square 
root of the products of bandpass amplitudes after fringe fit- 
ting. Examples of residual phases before and after calibration 
are shown in figure H) 

3.3. Computation of total group delays and phase delay 

rates 

The results of the fringe search are residual phase and group 
delay as well as their time derivatives determined from analy- 
sis of an observation at a given baseline of a given scan with 
respect to the a priori delay model used by the correlator. We 
need to compute the total path delay related to a certain mo- 
ment of time common to all observations of a scan, called a 
scan reference time (fsn)- 

For logistical reasons, fringe searching at different base- 
lines is performed independently and therefore, each obser- 
vation has its own reference epoch, called a fringe reference 
times (ffit). This time epoch is computed as the weighted mean 
epoch of a given observation. Since stations usually start and 
end at slightly different time and the number of processed 
samples may be different, in general, ffn is different at dif- 
ferent baselines of the same scan. The scan duration may be 
significantly different at different stations either by a schedule 
design when antennas with different sensitivities participate in 
observations — this is often made in geodetic observations, or 
due to losses of some of the observing time at stations for vari- 
ous reasons. If we set fsn as an average of ffn over all baselines 
of a scan, it may happen that for some baselines the difference 
ha- kit may reach several hundred seconds. Setting a com- 
mon fsrt for as many baselines as possible is desirable since 
it allows computation of delay triangle misclosures and some 
other important statistics. On the other hand, the uncertainty 
of the group delay estimate is minimal at ffn. The uncertainty 
in group delay at fsn grows as 

Cr^&rt) = O-rCffrt) + 2C0V(T,7-)(tsrt-tfrt) + CT^fet " tfe)' , (15) 



which is undesirable. In our work, we set the tolerance limit 
for the growth of the the uncertainty due to the differences 
between fsn and ffrt: O.Ict or 5 ps, whichever is less. Set- 
ting this limit, we find for each observation the maximum 
allowed Ifsri-ffrtl by solving quadratic equation [TS] In the 
case where all observations of a scan have overlapping inter- 
vals of acceptable scan reference times, we set it to the value 
that minimizes 2 ^ Cov(r , f) (tjit - tf,-t) + cr^iUn - tfn)^ over all 
observations. In the case where there are observations that 
have intervals of acceptable f,,, which are not overlapping, the 
set of observations of a scan is split into several subsets with 
overlapping acceptable fsn and the optimization procedure is 
performed under each subset. Finally, t^n is rounded to the 
nearest integer second. 

The VLBA correlator shifts the time tag of data streams 
from each station to the moment of time when the wavefront 
reaches the center of the coordinate system. This operation fa- 
cilitates correlation and allows station-based processing. The 
a priori path delay is computed for this modified time tag. 
This shift of the time tag depends on the a priori parameters 
of the geometric models, and therefore the total path delay 
produced from such a modified quantity would depend on er- 
rors of the a priori model which would considerably compli- 
cate data analysis. Therefore, we have to undo this shift of 
the time tag for the reference station of a baseline for further 
processing. 

The correlator delay model for the VLBA hardware and 
software correlators is computed as a sum of a 5th degree 
polynomial fit to the geometric delays, the linear clock off- 
set, and the coarse atmospheric model delay over intervals of 
120 s length. We find the a priori delay of the baseline refer- 
ence station r^f related to the time tag at TAI of the wavefront 
arrival to its phase center from the implicit equation 

k=5 

= 2^ fl/ [fsn - to - [Ta " " T,i r„ - T,,^ j j , (16) 

which is solved by iterations. Here to is the TAI time tag of 
the polynomial start time, and are the clock model 
and the atmosphere contribution of the a priori path model. 
We set on the right hand side of expression [16] to zero for 
the first iteration. Three iterations are sufficient to reach the 
accuracy 0. 1 ps. Using the value r^f found at the last iteration, 
we compute the a priori path delay for the remote station of 
the baseline: 

rr = T.a'r{tsn-to-{r^-r^-r:rr^-rZ))'- (H) 

The a priori delay rate is computed using an expression sim- 
ilar to equation[T7J Finally, we compute the total path delay by 
extrapolating the residual delay to the scan reference time: 

'''tot = ■^ap-'^^r + "^res + T^esitsn- tfa) ■ (18) 

The delay produced this way is the difference between the 
interval of proper time measured by the clock of the remote 
station between events of arrival of the wave front to the ref- 
erence point of the remote antenna and clock synchronization 
in TAI and the interval of proper time measured by the clock 
of the reference station between events of arrival of the wave 
front to the reference point of the reference antenna and clock 
synchronization in TAI. 
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Figure 4. Residual fringe phases (radians) before (left) and after (right) applying bandpass calibration in experiment bpl25b. 



3.4. Astrometric Analysis: Delay Modeling 

Our computation of theoretical time delays in general fol- 
lows the approach presented in detail by .Sovers et al. (199^ 
with some refinements. The most significant ones are the fol- 
lowing. The advanced expre ssion for time delay derived by 
iKopeikin and Schafed (119991) in the framework of general rel- 
ativity was used. The displacements caused by the Earth's 
tides were computed using the numerical values of the gen- 
eralized Love numbers presented by Mathews (2001) follow- 
ing a rigorous algorithm described at iPetrov and Ma (.2003i) 
with truncation at a level of 0.05 mm. The displacements 
caused by ocean loading were computed by convolving the 
Greens' func tions wi t h ocea n tide models using t he NLOADF 
algorithm of lAgnewl ( [19971) . The GOTOO model la^ dHH) 
of diurnal and semi -diurnal ocean tides, th e NA099 model 
of oce an zonal tides (M atsum oto et al.ll2000l) . the equilibrium 
model (iPetrov and Mai2003i) of the pole tide, and the tide with 
period of 18.6 years were used. Atmospheric pressure loading 
was computed by convolving the Greens' functions with the 
output of the atmosp here NCEP Reanalysis numerical model 
(iKalnay et al.lll996l) . The algorithm of comp utations is de- 
scribed in full details in lPetrov and BoyI (120041) . The empirical 
model of harmonic variations in the Earth orientation param- 
eters he o_2 9 1 2 10 derived fr om VLB I observations ac- 
cording to the method proposed bv lPetrovl (12007 a) was used. 
The time series of UTl and polar motion from the Goddard 
operational VLBl solutions were used as a priori. 

The a priori path delays in the neutral atmosphere in di- 
rections towards observed sources were computed by nu- 
merical integration of differential equations of wave prop- 
agation through the heterogeneous media. The four- 
dimensional field of the refractivity index distribution was 
computed using the atmospheric pressure, air temperature 
and specific humidity taken from the output of the Mod- 
ern Era Ret rospective- Analysis f or Research and Applications 
(MERRA) ( ISchubert et alj|2008l) . That model presents the at- 
mospheric parameters at a grid 1 /2° x 2/3° x 6'' at 72 pres- 
sure levels. 

For considering the contribution of the ionosphere to the 
phase of the cross-correlation spectrum, notice that the elec- 



tromagnetic wave propagates in a plasma with phase velocity 



(19) 



Nye"- 



where A^,, — electron density, e — charge of an electron, 
— mass of an electron, e,, — permittivity of free space, lu — 
angular frequency of the wave and c — velocity of light in 
vacuum. Phase velocity in the ionosphere is faster than the 
velocity of light in vacuum. 

After integration along the ray path and expanding expres- 
sion [19] withholding only the term of the first order, we get 
the following expression for additional phase rotation caused 
by the ionosphere: 



A0 = — , 

U! 

where ut is the angular frequency and a is 

-) 
e~ 



Icnie eo 



Nydsi- 



Ny dso 



(20) 



(21) 



Here si and ^2 are the paths of wave propagation from the 
source to the first and second station of the baseline. If J A^,. ds 
is expressed in units of 1 • 10'^ electrons/;7z^ (so-called TEC 
units), then after having substituted values of constants, we 
get a = 5.308018 • lO'" sec"' times the difference in the TEC 
values at the two stations. 

Since the ionosphere contribution is frequency-dependent, 
it distorts the fringe-fitting result. Taking into account that the 
bandwidth of the recorded signal is small with respect to the 
observed frequency, we can linearize eqn. [20| near the refer- 
ence frequency (j) = -a/uJo + a(uJi-uj)/Ldl. Comparing it 
with expression[4| we see that the first frequency-independent 
term contributes to the phase delay and the second term, lin- 
ear with frequency, contributes to the group delay. The fine 
fringe search is equivalent to solving the LSQ for Tp and Tg 
using the following system of equations: 



' Available at |http : / / astrogeo . org/erm| 



TpUJo + Tg(uJk - l^o) = (f>i + — , 



(22) 



where index ; runs over frequencies and index k runs over 
accumulation periods. 
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A solution of the 2x2 system of normal equations that orig- 
inates from equations l22l can be easily obtained analytically. 
Gathering terms proportional to a, we express the contribution 
of the ionosphere to phase and group delay as t^""" = -a/ cj^ 
and t]™° = fl/w- 
frequencies: 



where ujp and ujg are effective ionospheric 




\ 



(uJi-OJo) 



(23) 




n n 



i 



where w, is the weight assigned to the fringe phase at the /th 
frequency channel. 

They have a clear physical meaning: if the wide -band signal 
was replaced by a quasi-monochromatic signal with a group 
or phase effective ionosphere frequency, then the contribution 
to group or phase delay of the wide-band signal would be the 
same as the contribution of the quasi-monochromatic signal 
at those effective frequencies. 

For computing the contribution of the ionosphere, we used 
the total electron contents (TEC) maps from analysis of linear 
combinations of GPS observables made at two frequencies, 
1.2276 and 1.57542 GHz. Using GPS-derived TEC maps for 
data reduction of astronomic observations, first suggested by 
|Rps (2000), has became a traditional approach for data pro- 
cessing. Analysis of continuous GPS observations from a 
global network comprising 100-300 stations makes it feasi- 
ble to derive an empirical model of the total electron contents 
over the span of observations using the data assimilation tech- 
nique. Such a model is routinely delivered by GPS data anal- 
ysis centers since 1998. For our analysis we used t he TEC 
model provided by the GPS analysis center CODE dSchaed 
[19981) . The model gives values of the TEC in zenith direction 
on a regular 3D grid with resolutions 5° x 5° x 2''. 

For the purpose of modeling, the ionosphere is considered 
as a thin spherical layer with constant height H above the 
Earth's surface. The typical value of H is 450 km. In or- 
der to compute the TEC from GPS maps, we need to know 
the coordinates of the point at which the ray pierces the iono- 
sphere — point / in figure |5] First, we find the distance from 
the station to the ionosphere piercing point D = \AJ\ by solv- 
ing triangle OAJ. Noticing that |(M| = i?® and = R^+H, 
we immediately get 



„ . cos/i 
p = arcsm jj- 



1 + - 



D = Rn 



/2^(l-sin(£ + /3))+(^ 



(24) 



where E is the elevation of the source above the horizon. 

Then the Cartesian coordinates of point J are f+Ds. Trans- 
forming them into polar coordinates, geocentric latitude and 




Figure 5. Ray traveling from the source s to antenna A pierces the top of the 
ionosphere in the point J and the bottom of the ionosphere in the point /. The 
ionosphere is considered as a thin layer. P is the point on the Earth's surface 
beneath the ionosphere piercing point J. 

longitude, we get arguments for interpolation in the 3D grid. 
We used the 3-dimensional B-spline interpolation by expand- 
ing the TEC field into the tensor products of basic splines 
of the 3rd degree. Interpolating the TEC model output, we 
get the TEC through the vertical path \JI„\. The slanted path 
is \JIo\/ cos/3. Therefore, we need to multiply the ver- 
tical TEC by 1 / cos /?(£), which maps the vertical path delay 
through the ionosphere into the slanted path delay. Here we 
neglect the ray path bending in the ionosphere. We also ne- 
glect Earth's ellipticity, since the Earth was considered spher- 
ical in the data assimilation procedure of the TEC model. 

Combining equations, we get the final expression for the 
contribution of the ionosphere to path delay: 



4^Veff 



TEC- 



1 



cos P(E) 



(25) 



where /eff is the effective cyclic frequency, and the plus sign is 
for the group delay, and the minus sign is for the phase delay. 

Computation of the theoretical path delay and its partial 
derivatives over model parameters is made with using soft- 
ware VTE0. 

3.5. Astmmetric Analysis: Parameter Estimation 

Astrometric analysis is made in several steps. First, each 
individual 24 hour session is processed independently. The 
parameter estimation model includes estimation of 1) clock 
functions presented as a sum of a 2nd degree polynomial and 
a linear spline over 60 minutes; 2) residual zenith atmosphere 
path delay for each station presented as a linear spline; 3) co- 
ordinates of all stations, except a reference station; and 4) co- 
ordinates of the target sources. The goal of the coarse solu- 
tion is to identify and suppress outliers. The main reasons for 

* Available at |http : //astrogeo ■ org/vtd| 
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outliers are a) errors in determining the global maximum of 
the fringe amplitude during fringe search and b) false detec- 
tions. Both errors decrease with increasing the SNR. Because 
of this, we initially run our solution by restricting to SNR > 
6, and then restore good detections with SNR in the range of 
[5,6]. 

After identifying outliers and removing them from our solu- 
tion, we apply estimates of the parameters of the a posteriori 
model to outliers, which allows us to predict the path delay 
with accuracy better than 500 ps, except for sources that had 
fewer than two detections. Then we re-run the fringe search 
for outliers and restrict the search window to ±1000 ps with 
respect to predicted delay. We also lower the SNR detec- 
tion limit to 4.8, since the number of independent samples 
in the restricted search window, and therefore, the probabil- 
ity of false detection at a given SNR is less. This procedure 
allows to restore from 50 to 80% of observations marked as 
outliers in the previous step. The weights of observables were 

computed as w = 1 / ^/aJ+r\F), where do is the formal un- 
certainty of group delay estimates and r{b) is the baseline- 
dependent reweighting parameter that was evaluated in a trial 
solution to make the ratio of the weighted sum of the squares 
of residuals to its mathematical expectation to be close to 
unity using the technique similar to that used for fine fringe 
search. 

Finally, we run a global VLBI solution that uses all avail- 
able observations to date, 7.56 million, from April 1980 
through August 2010 in a single LSQ run. The estimated pa- 
rameters are 

— global (over the entire data set): coordinates of 4924 
sources, including target objects in VGaPS campaign, 
positions and velocities of all stations, coefficients of 
B-spline expansion that model non-linear motion of 17 
stations, coefficients of harmonic site position varia- 
tions of 48 stations at four frequencies; annual, semi- 
annual, diurnal, semi-diurnal, and axis offsets for 67 
stations. 

— local (over each session): tilts of the local symmetric 
axis of the atmosphere (also known as "atmospheric 
azimuthal gradients") for all stations and their rates, 
station-dependent clock functions modeled by second 
order polynomials, baseline-dependent clock offsets, 
and the Earth orientation parameters. 

— segmented (over 20-60 minutes): coefficients of linear 
splines that model atmospheric path delays (20 minutes 
segment) and clock functions (60 minutes segment) for 
each station. The estimates of clock functions absorb 
uncalibrated instrumental delays in the data acquisition 
system. 

The rate of change for the atmospheric path delays and 
clock functions between adjacent segments was constrained 
to zero with weights reciprocal to 1.1 • lO"''* and 2 • lO"''*, re- 
spectively, in order to stabilize solutions. We apply no-net 
rotation constraints on the positions of 212 sources marked 
as "defining" in the ICRF catalogue (iMa et al.lll998h that re- 
quires the positions of these source in the new catalogue to 
have no rotation with respect to the position in the ICRF cat- 
alogue to preserve continuity with previous solutions. 

The global solution sets the orientation of the array with 
respect to an ensemble of ^ 5000 extragalactic remote radio 



sources. The orientation is defined by the series of Earth ori- 
entation parameters and parameters of the empirical model of 
site position variations over 30 years evaluated together with 
source coordinates. Common sources observed in VGaPS as 
atmosphere and amplitude calibrators provide a connection 
between the new catalogue and the old catalogue of compact 
sources. 

3.6. Astrometric analysis: assessment of weights of 
observations 

As follows from the Gauss-Markov theorem, the estimate 
of parameters has minimum dispersion when observation 
weights are chosen reciprocal to the variance of errors. The 
group delays used in the analysis have errors due to the ther- 
mal noise in fringe phases and due to mismodeling theoretical 
path delay in the atmosphere: 

^' = 4 + 4 + '^L (26) 

where aJi^ is thermal noise, ct?, and a^^^ are the contribution of 
the ionosphere and the neutral atmosphere to the error budget. 

3.6.1. A priori errors of the GPS ionosphere model 

The first term, was estimated during the fringe fit- 
ting. The second term can only be evaluated indirectly. 
ISekido et al.l (l2003h used six dual-band intercontinental VLBI 
sessions at 10 stations in July 2000 and compared TEC val- 
ues from GPS with TEC estimated from VLBI observables. 
They drew a conclusion that the errors in path delay derived 
from the GPS TEC model were at the range of 70 ps at the 
zenith direction at 8.6 GHz. However, the ionosphere path 
delay is a non-stationary process. Therefore, great caution 
should be taken in an attempt to generalize conclusions made 
from analysis of a small network over a short time period. 
The non-stationarity of ionospheric fluctuations implies that 
there does not exist an exact expression for the variance of 
the ionosphere fluctuations during any given period, and any 
expression for the variance is an approximation. 

Since June 1998 when the GPS TEC maps became avail- 
able through August 2010, more than 2000 dual -band S/X 
VLBI sessions under geodesy and absolute astrometry pro- 
grams were carried out, including 92 sessions at the VLBA. It 
can be easily shown that the contribution of the ionosphere at 
X-band, Tv, can be found from the linear combination of group 
delay observables with coefficients that are expressed through 
effective ionosphere frequencies at these bands and ujs de- 
fined in expressionl23l 

Txi = {Tx - T,) , • , . (27) 

We used this dataset for evaluation of the errors of the 
contribution of the ionosphere to group delays derived from 
GPS TEC maps, considering the ionosphere contribution from 
dual-band VLBI observations as true for the purpose of this 
comparison. We computed the ionosphere contribution from 
the GPS model and from VLBI observations for each session. 
The root mean squares (rms) of the differences of the con- 
tribution VLBI-GPS was computed for all sessions and all 
baselines. We sought regressors that can predict rms(VLBI- 
GPS). We expect the short-term variability of the ionosphere 
at scales less than several hours to dominate the errors of the 
GPS model. The sparseness of the GPS network and limited 
sky coverage result in missing high frequency spatial and tem- 
poral variations of the ionosphere. The turbulent nature of the 
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Figure 6. Dependence of the rms of the differences VLBI minus GPS as a function of RG=^mis^j +rms^2 for baseline FD- VLBA/PIETOWN, 565 km long (left) 

and for baseline LA-VLBA/MK-VLBA, 4970 km long (right). The quantity rmSg, is the rms of GPS path delay at the i station of a baseline during a session. Each 
green dot corresponds to one VLBI session. The thick blue straight light is a linear fit through the data. 
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Figure 7. The floor of the regression model of the dependence of nTis(VLBI- 
GPS) of RG = ^rms^j +rms^, for all VLBA baselines as a function of the 
baseline length. 
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Figure 8. The slope of the regression model of dependence of rms(VLBI- 
GPS) of RG = 



baseline length. The straight lines show a linear approximation of the slope 
for two ranges of the baseline length: below and over 2000 km. 



ionosphere path delay variations suggests that the rms of the 
ionosphere model errors due to missed high frequency varia- 
tions will be related to the rms of the low frequency variations 
either as a linear function or as a power law. After many tries 
we found that the following parameter can serve as a regres- 

sor: RG = J rms^j +rms^2' where rmSg, is the rms of the GPS 



path delay at the / th station of a baseline during a session. The 
rmSgi and rmSg2 are highly correlated at short baselines and 

^rms^j +rms^2 > rms(g2 - gl). At long baselines the iono- 
sphere contribution de-correlates. Therefore, the dependence 
of rms(VLBl-GPS) versus RG will depend on the baseline 
length and possibly on other parameters. Figure |6] shows ex- 
amples of this dependence for a short baseline and for a long 
baseline. 

It is remarkable that rms(VLBl-GPS) versus RG fits rea- 
sonably well with a linear function. We computed rms(VLBl- 
GPS) for each baseline and fitted it to the linear model 
F + S ■ RG. The floor of the linear fit has a mean value around 
20 ps (Figure IT). As expected, the slope of the fit increases 
with baseline length as shown in Figure [8] The growth is 
linear up to baselines lengths of 2 000 km, which apparently 
correspond to the decorrelation of the paths through the iono- 
sphere. The growth of the slope beyond baseline lengths of 
2 000 km is slower and it shows more scatter. 

We use this dependence for predicting the rms of the GPS 
ionosphere model errors. For each station that participated in 
the experiment we computed the rms of the ionosphere varia- 
tions in zenith direction. Then we express the predicted vari- 
ance of the GPS ionosphere model errors as 



J I 



feff 



cos ^(£2) 



'""O'l cos^(£i) 
2 



(28) 



where r,„„o , is the ionosphere path delay at the ; th station 
computed using the GPS TEC maps, f \ the zenith ionosphere 
path delay from GPS TEC maps averaged over a 24* period 
with respect to the central date of the session, Fb and St are 
parameters of the linear model rms(VLBl-GPS) versus RG 
for a given baseline, /, is the frequency for which the model 
was computed (8.6 GHz), and /eft the frequency of the exper- 
iment for which the model is applied (24.5 GHz). The term 

Tiono.i — cosp(E) difference between the ionosphere path 

delay from GPS at a given direction and the average iono- 
sphere path delay scaled to take into account the elevation 
dependence. This term is an approximation of rms^i used for 
computation of RG. 

3.6.2. A priori errors of the path delay in the neutral atmosphere 
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Rigorous analysis of the errors of modeling the path delay 
in the neutral atmosphere is beyond the scope of this paper. 
Assuming the dominant errors of the a priori model are due to 
high frequency fluctuations of water vapor at scales less than 
3-5 hours, we seek the regression model in the form of de- 
pendence of the rms of errors on the total path delay in the 
non-hydrostatic component of the path delay. We made sev- 
eral trial runs using all 123 observing sessions at the VLBA 
under geodesy and absolute astrometry programs with recip- 
rocal weights modified according to 

'^L/ = ^'+(«-^) ■ (29) 

Here Ts is the contribution of the non-hydrostatic con- 
stituent of the slanted path delay and t~ is the non-hydrostatic 
path delay in the zenith direction computed by direct inte- 
gration of the equations of wave propagation through the at- 
mosphere using the refractivity computed from the MERRA 
model, and a is the coefficient. We found that when coeffi- 
cient a = 0.02 is used, the baseline length repeatability, de- 
fined as the rms of the deviation of baseline length with re- 
spect to the linear time evolution, reaches the minimum. We 
adopted value 0.02 in our analysis of VGaPS experiments. 
For typical values of t-, the added noise is 8 ps in zenith di- 
rection and 80 ps at 10° elevation. 

3.6.3. Ad hoc added variance of the noise 

We also computed for each baseline and each session an ad 
hoc variance of observables that, added in quadrature, makes 
the ratio of the weighted sum of squares of post-fit residuals to 
their mathematical expectation close to unity in a similar way 
as we updated fringe phase weights. Expression[T4lwas used 
for computing this variance. This ad hoc variance was applied 
to further inflate the a priori observable uncertainties that have 
already been corrected for the inaccuracy of the a priori model 
of wave propagation through the ionosphere and the neutral 
atmosphere according to expressions[28]and[29l In contrast to 
|28]and[29] the baseline-dependent ad hoc variance is elevation 
independent. 

4. VALIDATION OF THE WIDE-BAND FRINGE 
SEARCH ALGORITHM 

Using VIMA , we detected 327 target sources versus 136 
targets detected with AIPS, since the AIPS detection limit is 
lowered by a factor of \/8 = 2.83. Thus, the yield of the ex- 
periments was improved by the factor of 2.4 - a very signifi- 
cant improvement that well justified our efforts to create a new 
software package for processing astrometric observations. 

Another difference is that the results of processing the data 
with AIPS have ambiguities in group delay that are recipro- 
cal to the minimum difference between intermediate frequen- 
cies. This means that the group delays are r^ + ZTr,, where 
Ts = 76.923 ns and K is an arbitrary integer number At the 
same time, the results of processing with the wide-band fringe 
fitting algorithm do not suffer this problem. The reason for 
group delay ambiguities is that the narrow-band fringe fitting 
algorithm implemented in AIPS first coherently averages the 
data within each IF, and in the second step of fringe fitting 
it processes a rail of narrow-band signals. The Fourier trans- 
form that describes the dependency of the amplitude of the 
coherent sum of the cross-spectrum on group delay has a pe- 
riodicity that is reciprocal to the minimum frequency sepa- 
ration of IFs, 76.923 ns in our case. The wide-band fringe 




Group delay (ns) 

Figure 9. Coherent sum of cross-con'elation amplitude, normalized to its 
value at the maximum as a function of group delay shift with respect to the 
maximum for VGaPS observing sessions. The thick green line shows result 
of processing with 'PIMA. . The thin blue hne shows results of processing 
with AIPS. The pattern of the thin blue line repeats with a period of 76.923 ns. 



fitting algorithm does not average the spectrum. Therefore, 
the periodicity of the Fourier transform of the coherent sum 
of the cross-spectrum in the wide-band algorithm is equal to 
the spectral resolution, i.e. 2 mks for VGaPS experiments. 
Figure |9] illustrates this difference. 

The lack of group delay ambiguities has a profound effect 
on determining source positions with poorly known a priori 
coordinates. In the presence of group delay ambiguities, we 
had to solve for source positions first using less precise so- 
called narrow-band group delays determined by arithmetic 
averaging group delays computed for each IF independently. 
The accuracy of these source position estimates is often insuf- 
ficient to reliably resolve group delay ambiguities, especially 
in the presence of narrow-band group delay outliers. When 
the number of used observations is large, say more than 10, 
the data redundancy allows us to detect the presence of in- 
correctly resolved ambiguities and fix the problem. But if 
the number of observations is small, chances are the error in 
group delay ambiguity resolution will not be noticed. In the 
past, source position estimates made with less than 8 obser- 
vations were considered unreliable. The use of the wide-band 
fringe fitting algorithm eliminates this problem entirely. Re- 
processing the old data revealed that group delay ambiguities 
for a considerable number of sources were indeed resolved 
incorrectly, which resulted in source position errors as large 
as 4'! In contrast, the wide-band fringe fitting algorithm pro- 
vides reliable estimates of source positions using a minimum 
redundancy of 3 observations. 

Since the wide-band fringe search algorithm is new, we 
would like to be sure that the new algorithm does not intro- 
duce new systematic errors with respect to the old one. As 
a validation test, we re-processed a set of 33 VLBA abso- 
lute astrometry/geo desy experiments under the RDV program 
jPetrov et al.ll2009l) and 12 VLB A absolute astrom etry exper- 
iments under the K/Q program dLanvi et alj|2010l) . Each test 
experiment had a duration of 24 hours. 

The RDV experiments were observed on a global network, 
including all 10 VLBA stations, with dual-band receivers at 
8.6 GHz (X band) and 2.3 GHz (S band), with 4 IFs aflocated 
at X band and 4 IFs allocated at S band. Fringe fitting, out- 
lier elimination, re-weighting, and in the case of AIPS, group 
delay ambiguity resolution, were made independently using 
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Table 5 

Solution statistics from 33 global RDV sessions processed with AIPS and 
TKLMA . The statistics in the central column were computed using all 
observations. The statistics in the right column were computed using 
observations with SNR > 10. 





AIPS 


VIMA 
SNR„„„ = 5.0 


■piMA 
SNR„i„ = 10.0 


No. obs used 


467 769 


531 299 


472 717 


fit wrms 


18.40 ps 


21.22 ps 


19.65 ps 


No. sources 


776 


800 


773 


wrms A^coseo 


0.10 mas 


0.12 mas 


0.14 mas 


wrms Ae 


0.10 mas 


0.10 mas 


0.12 mas 


Bas. rep. at 5000 km 


4.81 mm 


4.75 mm 


4.96 mm 


Bas. rep. at 9000 km 


8.54 mm 


8.08 mm 


7.95 mm 



VIAiA and AIPS. Subsequent data reduction and parameter 
estimation were made using identical setups. Estimated pa- 
rameters of the solution were the same as in processing the 
VGaPS sessions, except for treatment of site positions: they 
were treated as local parameters, i.e. estimated for each ses- 
sion independently. 

The statistics of the solution for 33 global RDV sessions 
using AIPS and VIMA are shown in table |5] The weighted 
root mean squares (wrms) of the postfit residuals is larger in 
the VIMA solution for two reasons. First, the VIAiA solu- 
tion contain 14% more points, mainly with SNR's in the range 
[5.0, 10.0], that were undetected by the traditional AIPS al- 
gorithm. The errors of these observables are systematically 
larger We rerun the VIAiA solution and excluded all points 
with SNR either at X- or S-band less than 10. The difference 
in wrms postfit residuals was significantly reduced. The sec- 
ond reason is that the group delay formal errors were inflated 
in VIMA processing to make the ratio of postfit residuals of 
fringe fitting to its mathematical expectation close to 1 . This 
was not done for the AIPS solutions. 

Another test of goodness of the solution and possible pres- 
ence of system atic errors is the so -called baseline length re- 
peatability test (iPetrov et al.l 120091) . We computed the wrms 
of baseline length estimates with respect to the fitted linear 
model of their evolution with time. The dependence of the 
baseline length wrms with the length of baselines L is fitted by 
function R{L) = VA^ + B^L. Values of R(L) at L=5000 km and 
L=9000 km are presented in table |5] We also computed the 
wrms of the deviations of estimates of daily offsets of nutation 
angles in longitude, AV-', and nutation in obliquity, Ae, with 
respect to the empirical nutation expansion he o_2 00 912 01. 
The statistics in table |5] show the satisfactory agreement be- 
tween AIPS and VIMA solutions and do not reveal any sys- 
tematic differences. 

Since the goal of VGaPS was to derive source positions, 
comparison of the positions from AIPS and VIMA process- 
ing is important to evaluate the level of systematic differences. 
We, thus, computed the differences in sources coordinates 
AacosS(a) and A6{6). We restricted our analysis to 528 
sources that had more than 64 observations in both AIPS and 
VIMA in order to avoid effects of a greater number of ob- 
servations available in the VIMA solutions. Plots of these 
differences are shown in figure [TO] and position comparisons 
are shown in table |6] The table contains the results from 528 
sources observed in RDV sessions with at least 64 used ob- 
servations. The 2nd column shows statistics of differences 
between AIPS and VIMA with an SNR cutoff of 5. The 3rd 
column shows the differences between AIPS and VIMA so- 



Table 6 

Differences between AIPS and VIMA positions of 528 sources observed 

in RDV experiments. The differences in the left column were computed 
using all observations. The differences in the right column were computed 
using observations with SNR > 10. 





AlPS-VIMA 


AWS-VIMA 




SNR„,„ = 5.0 


SNR„„„ = 10.0 


wrms Aacos(5(Q) 


0.072 mas 


0.020 mas 


wrms A(5((5) 


0.088 mas 


0.032 mas 



lutions with an SNR cutoff of 10. These differences are only 
about 20% of a typical position uncertainty. 

In the right plot of Figure [TO] we can see a small system- 
atic declination differences at low declinations. Similar differ- 
ences but produced from VIMA solution made with the SNR 
cutoff 10.0, shown in Figure [TT] help to understand the origin 
of this pattern. When the SNR cutoff is raised to 10, the wrms 
of differences are reduced by 2-3 times to 0.02-0.03 mas and 
the systematic pattern disappear. For comparison, the aver- 
age formal uncertainties for declinations are 0.12 mas and for 
right ascensions scaled by cos (5 are 0.07 mas, so the system- 
atic errors are small. Sources at low declinations observed on 
a VLBl array located in the northern hemisphere are neces- 
sarily taken at low elevations. Including in a solution addi- 
tional low SNR observations at low elevations, unavailable in 
the AIPS solutions, changes the contribution of systematic er- 
rors due to mismodeling path delay in the neutral atmosphere 
that are higher at low elevations. At present, it is not clear 
whether including observations with SNR in range 5-10 in- 
creases systematic errors, or the opposite, including these ob- 
servations reduces the systematic errors. However, the mag- 
nitude of the differences, less than 0.5 mas as declinations in 
the range [-50° -25°] does not raise a concern. 

Finally, we have re-analyzed 12 VLBA experiments under 
the K/Q program. The frequency setup in K/Q and VGaPS 
campaigns was identical. Lowering the detection limit by a 
factor of \/8 w 2.83 with the use of the wide-field fringe fitting 
algorithm greatly helped. Processing the data with VIMA 
did not detect only 8 out of the 340 observed sources, with 61 
non-detections in data processing with AIPS. We compiled 
table |7] of statistics of the K/Q solution in a form similar to 
table |5] Analysis of source position differences did not reveal 
any pattern of systematic errors. The statistics of these differ- 
ences presented in table [8] do not exceed formal uncertainties 
of positions which are 0.08 mas and 0.14 mas in right ascen- 
sions scaled by cos 5 and declination, respectively. The first 
column is four sources with an SNR cutoff of 5, the second 
column with an SNR cutoff of 5 • VS = 14.1. The baseline 
length repeatability and the wrms of nutation offset time se- 
ries from VIMA solutions are 20-30% smaller We do not 
have an explanation why VIMA solution produces notice- 
ably better results for 24 GHz observations, but no significant 
improvement was found from analysis of S/X RDV experi- 
ments. 

The above results of analysis of validation runs processing 
0.6 million observations at K, X and S bands, collected during 
1080 hours of recording at the VLBA and the global VLBl 
network with both AIPS and VIMA demonstrate that the 
new wide-band algorithm for group delays does not introduce 
any significant systematic errors while detecting more sources 
because it evaluates group delays using the coherent sum of 
the data across the wide-band. We conclude that VIMA has 
passed major vaUdation tests. 
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Figure 10. Systematic differences in source coordinate estimates between solutions using group delays derived by 'PIA4A and by AIPS. Left plot shows 
A« cos (5 (a), right plots shows A5(5). 

5. INVESTIGATION OF SYSTEMATIC ERRORS IN 
SOURCE POSITIONS 

Using single-band data, ionosphere path delay mismodeling 
may produce systematic position errors. iLanvi et al] (120101) 
showed that in their analysis of K-band observations system- 
atic errors reached several mas and ha d a tendency to be larger 
at low declinations. In section [J.6.1l we evaluated the rms of 
random errors caused by ionosphere path delay mismodeling. 
However, inflating weights to account for the variance of er- 
rors in general does not guarantee that source positions will 
have no systematic errors. 

To evaluate the magnitude of possible ionosphere driven 
systematic errors we made the following Monte Carlo sim- 
ulation. We added to the theoretical path delay the zero-mean 
Gaussian noise M-A^(0,(7,), with the variance ct/ computed ac- 
cording to expressionl28] The noise was magnified M times in 
order to make the contribution of ionosphere path delay errors 
dominant over other sources of errors. We made 64 analysis 
runs of VGaPS data using different seeds of the random noise 
generator The magnification factor 100 was used. Thus, we 
produced 64 estimates of position of each source with added 
noise. We computed the rms of position estimates of each tar- 
get source and divided it by M. To the extent of the validity 
of expressionism these rms's represent expected errors due to 
the inadequacy of the ionosphere path delay models based on 
using the TEC models derived from GPS analysis. Plots of 
Aq!,((5) and A(5,((5) errors are shown in Figures fTSHTSl For 
90% of the sources, errors are at the level of 0.02-0.04 mas. 
A6(6) increases to 0.1 mas at declinations less than -20°, and 
for some sources may reach 0.4 mas. The disparity in system- 
atic errors for sources at comparable declinations reflects the 
disparity in the number of observables used in the solution. 

For assessment of remaining systematic errors we exploited 
that fact that 56 known sources were observed as amplitude 
and atmospheric calibrators. Positions of these sources are 
known from previous dual-band S/X observations with accu- 
racies better than 0.1 mas. We split the set of 56 calibrators 
into two subsets of 28 objects and ran two additional solu- 
tions. In the first solution we suppressed 28 calibrators in all 
sessions except VGaPS and determined their positions solely 
from VGaPS. In the second solution we did the same with the 
second subset. Considering that the positions of calibrators 
from numerous S/X observations can be regarded as true for 
the purpose of this comparison, we treated the differences as 
VGaPS errors. 
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Figure 11. Systematic differences A5(<5) in source coordinate estimates be- 
tween solutions using group delays derived by 'PIMA, with SNR cutoff 10 
and by MPS. 



Table 7 

Solution statistics from 12 VLBA experiments at 24 GHz under the K/Q 
program processed with AIPS and 'PXMA . The statistics in the central 
column were computed using all observations. The statistics in the right 
column were computed using observations with SNR > 1 0. 





AIPS 


PIMA PIMA 
SNR„„„ = 5.0 SNR„i„=14.1 


# Obs used 


104 887 


139213 


105 761 


fit wrms 


19.50 ps 


23.08 ps 


19.19 ps 


# Sources 


279 


332 


276 


wrms A^coseo 


0.13 mas 


0.10 mas 


0.09 mas 


wrms Ae 


0.18 mas 


0.13 mas 


0.14 mas 


Bas. rep. at 5000 km 


5.56 mm 


4.58 mm 


4.58 mm 


Bas. rep. at 9000 km 


9.29 mm 


7.36 mm 


7.32 mm 



Table 8 

Differences between AIPS and PIMA positions of 244 sources observed 

in the K/Q VLBA experiments. The differences in the left column were 
computed using all observations. The differences in the right column were 
computed using observations with SNR > 10. 





AWS-PIMA 


AWS-PIMA 




SNR,„„ = 5.0 


SNR„i„ = 10.0 


wrms Aa cos (5(a) 


0.062 mas 


0.060 mas 


wrms A<5((5) 


0.108 mas 


0.096 mas 
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Figure 12. Modeled systematic errors Aa, (5) cos <5 driven by tiie mismodel- 
ing ionospiiere patli delay contribution evaluated from the Monte Carlo sim- 
ulation. 
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Figure 13. Modeled systematic eiTors A5'(5) driven by the mismodeling 
ionosphere path delay contribution evaluated from the Monte Carlo simula- 
tion. 



We computed the per degree of freedom statistics for the 
differences in right ascensions and declinations Aa and AS 
and sought additional variances and vg that, being added 
in quadrature to the source position uncertainties, will make 
them close to unity: 

k=n 



Aal cos^ 6k 



lidF 



k=n 



k=l 



k=n 



(30) 



ndF 



k=i 



k=n 



k=l 



The denominator in equation [30] is a mathematical expec- 
tation of the sum of squares of differences, provided the esti- 
mates of source positions are statistically independent. 

We found the following additive corrections of the uncer- 
tainties in right ascensions scaled by cos 5 and for declinations 
respectively: = 0.08 mas and V5=0.120 mas. 

The final inflated errors of source positions, al^if) and 



crl{f) = <Tl + vl + aj/cos^S 



(31) 



Position of so me sources may indeed be different d ue to the 
core-shift effect (iKovalev et alJ 120081 lPorcasll2009h . Treat- 
ment of this shift as VGaPS errors makes our estimates of 
reweighting parameters and therefore reported final inflated 
errors somewhat too conservative. 

6. IMAGING 

For imaging purposes we performed the a priori calibration 
of the data following a traditional way, using MPS do reisen I 
119881) . In the future, we plan to introduce all extra steps re- 
quired for an accurate amplitude calibration into TZLMA as 
well. 

We followed the usual AIPS initial VLBA calibration pro- 
cedure involving a priori amplitude calibration with mea- 
sured antenna gain curves and system temperatures as well 
as sampling-based calibration adjustments. Atmospheric ab- 
sorption is significant at 24 GHz. We have estimated its effect 
using system temperature data covering the whole range of 
elevations and weather information in order to adjust visibil- 
ity amplitudes for the opacity. Typical values of the opacity 
were found to range between 0.03 and 0. 1 for different VLBA 
telescopes and observing epochs. We performed phase cali- 
bration using the phase calibration signal injected during ob- 
servations and fringe fitting. A separate solution for station- 
based group delay and phase delay rate was made for each 
frequency channel (IF). As the final step of calibration, band- 
pass corrections were determined and applied. 

Our observations were scheduled around 24 GHz since the 
K-band continuum performance is better at this frequency, 
away from the water line. Unfortunately, at the time of these 
observations most of the VLBA telescopes had no gain curve 
measurements close to 24 GHz. It has changed since 2007 
when regular 24 GHz gain curve measurements started to be 
performed at all VLBA stations. For all VLBA antennas but 
MK-VLBA and HN-VLBA we have used gain curves measured 
at 22.2 GHz while for the former the curves at 23.8 GHz were 
applied. Antenna efficiency as well as t he noise diode spec- 
trum do change with frequency (see, e.g.. lPetrov et alj|2007d 
as well as results of VLBA gain curve measurements at 22 
and 24 GHz after 200lQ). This is one of the main sources of 
the total amplitude calibration uncertainty. Additionally, IFs 
in our experiment are widely spread (Table[2]i which might in- 
troduce extra amplitude shifts. We used strong flat-spectrum 
sources in the sample in order to estimate global relative am- 
plitude correction factors for different IFs but did not find any 
to be larger than 10% with high confidence. No extra fre- 
quency channel specific amplitude corrections were applied 
to the data. 

After a priori calibratio n, data were i mported to the 
Caltech DIFMAP package (IShepherdlll997b . visibility data 
flagged, and images were produced using an automated hy- 
brid inmging_procedure originally suggested by Greg Tay- 
lor ( iPearson et al.lll99"4h which we optimized for our exper- 
iment. The procedure performs iterations of phase and ampli- 
tude self-calibration followed by CLEAN image reconstruc- 
tion. We were able to reach the VLBA image thermal noise 
level for most of our final CLEAN images. Examples of 3 

"^ [http : / /www . vlba . nrao . edu/ astro/VOBS/ astronomy/vlba_gains . kj 



18 



Petrov et al. 



clev=0.0012, max = 0.145 Jy/beom. freq = 24 GHz 



clev=0.0014, mQx = 0.349 Jy/beam. freq = 24 GHz 



clev = 0.0018, max = 0.239 Jy/beam, freq = 24 GHz 
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Figure 14. From top to bottom. Row 1: Naturally weighted CLEAN images at 24 GHz. The lowest contour is plotted at the level given by "clev" in each panel 
title (Jy/beam), the peak brightness is given by "max" (Jy/beam). The contour levels increase by factors of two. The dashed contours indicate negative flux. The 
beam is shown in the bottom left corner of the images. Row 2: Dependence of the correlated flux density on projected spacing. Each point represents a coherent 
average over one 2 min observation on an individual interferometer baseline. The error bars represent only the statistical errors. 



images, for compact and resolved objects, are shown in Fig- 
ure 111 

Total errors of our measurements of correlated flux density 
values for sources stronger than ^--^200 mJy were dominated 
by the accuracy of amplitude calibration described above. 
This considers the error of amplitude calibration as not ex- 
ceeding 15 % and this estimate is confirmed by our compar- 
ison of the flux densities integrated over the VLBA images 
with the single-dish flux densities which we measured with 
RATAN-600 in 2006 for slowly varying sources without ex- 
tended structure. Details of the RATAN-600 single-dish ob- 
serving program including the method of observat i ons a nd 
data processing can be found in lKovalev et"aLl (119991 120021) . 

7. THE CATALOGUE OF SOURCE POSITIONS 

Of 487 sources observed, three or more detections were 
found for 327 objects. After careful identification and re- 
moval outliers due to the incorrect selection of the global max- 
imum for weak sources with SNR < 6, we selected 33,452 
observations out of 59,690 from three VGaPS experiments for 
analyzing in the single LSQ solution together with 7.56 mil- 
lion other VLBI observations. The semi-major error ellipses 
of inflated position errors for all sources except 0903-1-154 
vary in the range 0.21 to 20 mas with the median value of 
0.85 mas. The histogram of position errors is shown in Fig- 
ure [B] 

The VGaPS catalogue is listed in Table |9] Although posi- 
tions of all 5047 astrometric sources were adjusted in the LSQ 
solution that included the VGaPS sources, only coordinates of 
327 target sources observed in the VGaPS campaign are pre- 



sented in the table. The first column gives calibrator class; 
"C" if the source is recommended as a calibrator or "U" if it 
has an unreliable position, since there were less than 5 detec- 
tions and there is a risk that the secondary maximum of the 
coherent sum of weighted complex cross-correlation samples 
has been picked and has not been flagged out. The second and 
third columns give the IVS source name (B 1950 notation) and 
lAU name (J2000 notation). The fourth and fifth columns give 
source coordinates at the J2000.0 epoch. Columns /6/ and /7/ 
give source position uncertainties in right ascension and decli- 




3 4 5 

Position uncertainty (mas) 

Figure 15. The histogram of the semi-major axes of inflated position error 
ellipses among 327 target sources in the VGaPS catalogue. The last bin shows 
errors exceeding 4.75 mas. 
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Table 9 

The VGaPS catalogue 



Source name J2000.0 Coordinates Errors (mas) Con'elated flux 

density (in Jy) 



Class 


IVS 


lAU 


Right ascension 


Declination 


Aa 


AS 


Corr 


#Obs 


Total 


Unres 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


C 


2358+605 


JOOOl+6051 


00 01 07.099852 


+60 51 22.79800 


0.94 


0.51 


-0.115 


123 


O.ll 


0.11 


C 


2359+606 


J0002+6058 


00 02 06.696680 


+60 58 29.83950 


4.42 


1.92 


-0.168 


26 


0.09 


< 0.08 


c 


0002+541 


J0005+5428 


00 05 04.363368 


+54 28 24.92414 


0.60 


0.47 


0.167 


79 


0.34 


0.11 


c 


0003+505 


J0006+5050 


00 06 08.249784 


+50 50 04.41150 


0.68 


0.81 


-0.041 


64 


0.16 


0.14 


c 


0005+568 


J0007+5706 


00 07 48.468649 


+57 06 10.43705 


1.81 


2.18 


0.631 


46 


0.08 


0.09 


c 


0012+610 


J0014+6117 


00 14 48.792125 


+61 17 43.54198 


0.57 


0.29 


-0.034 


153 


0.25 


0.16 


c 


0024+597 


J0027+5958 


00 27 03.286191 


+59 58 52.95899 


0.56 


0.34 


-0.214 


136 


0.23 


0.16 


c 


0032+612 


JO035+613O 


00 35 25.310605 


+61 30 30.76122 


0.94 


0.52 


-0.038 


117 


0.13 


0.10 


c 


0034-220 


J0037-2145 


00 37 14.825799 


-21 45 24.71171 


1.17 


2.78 


-0.834 


59 


0.09 


0.09 


c 


0039+568 


J0042+5708 


00 42 19.451680 


+57 08 36.58569 


0.39 


0.25 


0.046 


162 


0.48 


0.32 


c 


0041+677 


J0044+6803 


00 44 50.759596 


+68 03 02.68540 


0.67 


0.29 


-0.163 


154 


0.23 


0.19 


c 


0044+566 


J0047+5657 


00 47 00.428864 


+56 57 42.39373 


0.53 


0.39 


0.006 


154 


0.18 


0.13 



Note. — Tablel9]is presented in its entirety in the electronic edition of the Astronomical Journal. A portion is shown here for guidance regarding its 
form and contents. Units of right ascension are hours, minutes and seconds. Units of dechnation ai"e degrees, minutes and seconds. 



nation in mas after applying inflated errors according to equa- 
tion |3T] (without cos 6 factor), and column /8/ gives the cor- 
relation coefficient between the errors in right ascension and 
declination. The number of group delays used for position de- 
termination is listed in column /9/. Column /lO/ gives the es- 
timate of the flux density integrated over the entire map. This 
estimate is computed as a sum of all CLEAN components if 
a CLEAN image was produced. If we did not have enough 
detections of the visibility function to produce a reliable im- 
age, the integrated flux density is estimated as the median of 
the correlated flux density measured at projected spacings less 
than 70 MA. The integrated flux density is the source total 
flux density with spatial frequencies less than 12 MA filtered 
out, or in other words, the flux density from all details of a 
source with size less than 20 mas. Column /1 1/ gives the flux 
density of unresolved components estimated as the median of 
correlated flux density values measured at projected spacings 
greater than 400 MA. For some sources no estimates of the 
unresolved flux density are presented, because either no data 
were collected at the baselines used in calculations, or these 
data were unreliable. 

The on-line version of this catalogue is posted on the Wet0- 
For each source it has 4 references: to a FITS file with 
CLEAN components of naturally weighted source images; 
to a FITS file with calibrated visibility data; to a postscript 
map of a source; and to a plot of correlated flux density as a 
function of the length of the basehne projection to the source 
plane. 

In Table [TOl we present a priori coordinates, total flux den- 
sities extrapolated to 24 GHz and spectral index estimates for 
160 target objects that have not been detected in the VGaPS 
experiment. Some of these sources were detected in other 
VLBl astrometry experiments at S/X bands. 

8. COMPARISON OF K- AND X/S-BAND 
ASTROMETRIC VLBl POSITIONS 

We searched the VLBl archive and found that among our 
target sources, 206 were observed with S/X at the VLBA un- 
der VCS and RDV programs before November 2010. We in- 
vestigated the differences in K-band observations against in- 
dependent S/X observations. We restricted our analysis to 192 
objects that had uncertainties from X/S and K band solutions 

^ |http : //astrogeo ■ org/vgaps| 



Table 10 

The list of 160 sources that have not been detected in VGaPS observations. 



Source names 


Right ascens. 


Declination g 


,al. lat. 


Flux Sp.Ind. 


# 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 






h m s 


o / // 


deg 


mJy 






2359+548 


J0002+5510 


00 02 00.470 


+55 10 38.00 


-6.8 


121. 


-0.03 


2 


0003+669 


J0006+6714 


00 06 10.000 


+67 14 38.30 


5.0 









0009+655 


J0012+6551 


00 12 37.671 


+65 51 10.82 


3.5 


195. 


-0.59 


8 


0010+722 


J0013+7231 


00 12 58.750 


+72 31 12.76 


10.1 


501. 


0.02 


8 


0017+590 


J0020+5917 


00 20 24.550 


+59 17 30.50 


-3.1 


324. 


-0.02 


6 


0018-194 


J0021-1910 


00 21 09.370 


-19 10 21.30 


-79.6 









0028+592 


J003 1+5929 


00 31 03.120 


+59 29 45.30 


-3.0 


3855. 


1.52 


2 


0041+660 


J0044+6618 


00 44 41.300 


+66 18 42.00 


3.7 


119. 


-0.91 


7 


0107+562 


JO 110+5632 


01 10 57.553 


+56 32 16.93 


-5.9 


267. 


-0.66 10 


0113+241 


JO 116+2422 


01 16 38.067 


+24 22 53.72 


-37.8 


168. 


-0.03 


7 


0128+554 


J0131+5545 


01 31 13.860 


+55 45 13.20 


-6.4 


150. 


-0.06 


5 



Note. — Table [To] is presented in its entirety in the electronic edition of the Astro- 
nomical Journal. A portion is shown here for guidance regarding its form and contents. 
Columns 1-2 show lAU B1950 and J2000 source names, columns 3-5 show a priori 
sources positions at the J2000.0 epoch, column 6 shows extrapolated a priori total flux 
density at 24 GHz, column 7 shows coarse estimate of the spectral index, and column 
8 shows the number of measurements of flux density found in the CATS database that 
were used for evaluation of the flux spectral index and extrapolation the flux density. 

less than 5 mas. Figures [16] show the differences in right as- 
censions and declinations. 

8.1. Special cases: 3C 1 19 and 3C410 

Positions of two sources, J0432H-4138 (also known as 
3C1 19) and J2020H-2942 (3C410) are found to be significantly 
off, namely 38.3 and 38.0 mas respectively. The significance 
of this offset, 55 and 72 times of inflated uncertainties, is too 
high to be explained by known error sources. RATAN-600 
observations have shown that both of them continuously show 
steep radio spectra and are slowly variable. On VLBl scales, 
they we re foun d to have significantly extended structure (e.g., 
Figures [T7TT8] l. 

One of these compact steep spectrum radio sources, the 
object J0432H-4138, is well studied at parsec scales. Re- 
cent high-dynamic range images at 5 and 8.4 GHz and 
full bibliography on his toric observations can be found in 
iMantovanni et al] (120101) . This source shows several bright 
components at 24 GHz. Two of them, namely components 
A' and 'C (Figure [TtIi. are located 40.6 mas apart. The fea- 
ture C is stronger, its total flux density is 0.54 Jy, but more 
extended. The FWHM of a circular Gaussian component fit- 
ted to the feature is found to be 1.4 mas. The feature A is 
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Figure 16. The differences in right ascensions scaled by cos 5 (top) and de- 
clinations (bottom) from K band VGaPS versus S/X historical VLBA obser- 
vations among 164 sources with uncertainties less than 2 mas. Positions of 
two sources differ significantly. 

dimmer, 0.13 Jy, but more compact — 0.3 mas. We note 
that here and below in this section the total flux density re- 
ported for different source features is calculated as a sum of 
all CLEAN components representing the corresponding struc- 
ture. The Gaussian components fit is being performed in the 
visibility plane. 

It is evident that the VGaPS 24 GHz observations referred 
the position of the weaker A component, while the S/X ob- 
servations referred the position of the C component. It is 
counter-intuitive. Since the component C is resolved, its con- 
tribution to fringe amplitude is small at long baseline projec- 
tions. At short baseline projections the component C dom- 
inates, at long baseline projections the component A domi- 
nates. But since the partial derivative of group delay with re- 
spect to source position is proportional to the baseline length 
projection on the source tangential plane, the contribution of 
long basel ine dominates in estimate of source position. Ac- 
cording to iMantovanni et alJ (1201 Oh . the total flux density of 
components A and C at 8.4 GHz in 2001 was 70 mJy and 
1121 mJy respectively. This large difference between compo- 
nents' flux densities is also supported by the previous VCSl 
observations of this object at S/X-bands. The component C 
is less resolved at 8.4 GHz and it dominates over A even at 
long baselines. Therefore, coordinate of a source at X-band 
are closer to the C component. We checked it directly, by sup- 
pressing observations at baselines longer than 800 km in a K- 
band trial solution. The position estimate became close to the 
X/S position. It is worth to note that the difference in position 
at K- and X-bands are 29.2 ± 0.4 mas in right ascension and 




Relative Right Ascension (nnas) 

Figure 17. Naturally weighted K-band VGaPS CLEAN image of 3C119. 
The lowest contour of 4. 1 mJy/beam is chosen at three times the rms noise, 
the peak brightness is 97 mJy/beam. The contour levels increase by factors of 
two. The dashed contours indicate negative brightness. The beam's full width 
at half maximum (FWHM) is shown in the bottom left corner of the images in 
grey. Red and blue spots indicate the positions and sizes (FWHM) of circular 
Gaussian model components for the features 'A' and 'C respectively. 



24.1 ±0.4 mas in declination, while the offset of the compo- 
nent A with respect to the component C at K-band is slightly 
larger: 30.25 ±0.10 mas and 27.05 ±0.15 mas. The K-band 
and X-band position is not exactly the position of components 
A and C, since in both solutions the contribution of the second 
component is small but not entirely negligible. 

The nature of the difference in positions of J2020-I-2942 is 
simila r. It was observed in VCS2 experiment on May 15, 
2002 ( iFomalont et al.ll2003l) . It had 63 detections at S-band 
and 72 detections at X-band. Position of this source at S-band 
from analysis of only S-band data with applying the iono- 
sphere contribution from the GPS TEC model shows a very 
large offset of 0."5 with respect to the X-band position (re- 
fer to Table fTTl i. The errors of the ionosphere contribution 
at S-band during the solar maximum affected position esti- 
mates considerably, however not to that extent. Comparison 
of positions of 130 sources with S > from the solution that 
used only S-band group delay observables with respect to the 
X/S solution in that experiment showed the differences in the 
range 2-7 mas. It is remarkable that the X/S position is away 
from both X and S-band positions, although intuitively we 
can expect them to be between X and S positions. This can be 
explained if to surmise that the source J2020-I-2942 has two 
components, 0".48 apart, one of them visible at X-band, but 
not visible in the original S-band image, and another visible 
at S-band, but not at X-band. An ionosphere-free linear com- 
bination of X and S band observables is used in the X/S solu- 
tion: ( 1 + /3) T^.r -fiTgs, where /3 = 1 / (w J /cj^ - 1 ) as it follows 

from eqn. |27] In the case if position at S-band, k^, is shifted 
with respect to the X-band position vector, k^, the ionosphere 
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Table 11 

Position difference of J2020+2942 at different frequencies with respect to its position at 

X-band 



Band 


RA shift 


DEC shift 


X 


0.0 ±0.58 mas 


0.0 ±0.71 mas 


S 


25.77 ±4. 11 mas 


478.79 ±4.39 mas 


X/S 


-2.29 ±0.54 mas 


-36.78 ±0.66 mas 


K 


0.25 ±0.41 mas 


1.45 ±0.70 mas 



Note. — The offset in right ascension is scaled by cos 5. 

free linear combination can be written as 

Vf = il + f3)Tg,,{l)- f3Tg,il) - p^il-l). (32) 

ok 

The first two terms correspond to the case of no offset be- 
tween X and S band positions. Therefore, estimates of the 
source position from the ionosphere-free linear combination 
of observables will be shifted at -/3 with respect to the off- 
set {ks-kx), i.e. in the opposite direction. Parameter (5 was 
0.076153 in the VCS2 experiment. Therefore, if our hypothe- 
sis that the X-band and S-band observations detected emission 
from two components is true, then the shift of X/S position 
with respect to the X-band position should be -1.96 mas in 
right ascension and -36.46 in declination, just within 0.3 mas 
from reported X/S positions! The K-band position is within 
1 .5(7 from the X-band position. 

In order to check our hypothesis, we have re-imaged VCS2 
observations of J2020H-2942 in a wide field and have detected 
several previously unknown features 'B', 'C, 'D', on a dis- 
tance up to about 500 mas from the dominating extended 
structure 'A' at S-band (Figure [TsTl. The total flux density of 
the features A and D is 1.40 Jy and 0.09 Jy respectively. The 
feature A is significantly more extended then D. We have fit- 
ted two circular Gaussian components to the Mv-data in order 
to determine positions of the features A and D. We note that 
the accuracy of position determination for the component A 
from the image is very poor since its structure being extended 
over at least 40 mas is not well represented by a single Gaus- 
sian component. The distance between Gaussian components 
A and D is 479 mas while positional difference in right ascen- 
sion and declination is 25 and 479 mas respectively — in a 
very good agreement with independent astrometric measure- 
ments (Table [TTTi. 

X-band and K-band wide field imaging (Figure [T9]l did not 
reveal components with wide separation. We conclude from 
the astrometric analysis presented above that X- and K-band 
images represent the more compact component D. In this case 
we could also analyze its spectrum on the basis of simultane- 
ous S/X-band observations. Its total flux density at X-band is 
found to be 0.25 Jy which provides the 2.3-8.6 GHz spectral 
index estimate a = +0.8 (flux density cx z^") — an indication 
of synchrotron emission with significant self-absorption. The 
features A, B, and C become too weak and/or too resolved 
that we could detect then in the snapshot VLBA images with 
a limited dynamic range and wv-coverage. 

8.2. General comparison: uncertainties, systematic 
K-S/X-band difference, and the core-shift effect 

Position differences for other objects do not show pecu- 
liarities. For instance, no declination de pendent sys t ematic 
differences similar to those reported by iLanvi et al] ( 120101) 
were found. The wrms of the differences is 0.46 mas in right 
ascension scaled by cos (J and 0.61 mas in decUnation. We 
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Figure 18. Naturally weighted S-band CLEAN image of 3C410 redone by 
us using VCS2 VLBA observations on 2002-05-14. The lowest contour of 
2.5 mjy/beam is chosen at three times the rms noise, the peak brightness is 
242 mJy/beam. The contour levels increase by factors of two. The dashed 
contours indicate negative brightness. The beam's FWHM is shown in the 
bottom left corner of the images in grey. Red and blue spots indicate the 
positions and sizes (FWHM) of circular Gaussian model components for the 
features 'A' and 'D' respectively. It should be noted that feature A is very 
extended and the single Gaussian component does not represent it well. 



have computed the normalized distances by dividing them by 
^Jej + el^, where e^ is the projection of the error ellipse of 
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Figure 19. Naturally weighted X-band (left, VCS2 data) and K-band (right, 
this survey — VGaPS) CLEAN images of 3C410. The lowest contour at 
X- and K-band is 3.3 and 2.8 mJy/beam while the peak brightness is 56 and 
178 mJy/beam respectively. On the basis of our analysi s, w e identify the 
feature presented in this Figure as feature 'D' from Figure [T8l 
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Figure 20. The empirical distribution of 190 normalized distances between 
K band position of target sources and their X/S positions (broken line) and 
the best fit Rayleigh distribution with cr = 0.90 (solid thick line). 



the K-band position to the direction of position difference and 
e„ is the similar projection of the error elHpse of the X/S po- 
sition. In the case if position errors from K-band and X/S 
catalogues are independent and Gaussian with the variance 
equal to reported uncertainties, the distribution of normalized 
distances will be Rayleigh with cr = 1 . The average of the nor- 
malized distances over all sources, except J0432-I-4138 and 
J2020H-2942, is 1.276, only 2% greater than the mean of the 
Rayleigh distribution, \pnjl. However, a close examination 
of the distribution (see Figure l20li reveals a slight deviation 
of its shape from the shape of the Rayleigh distribution. The 
Rayleigh distribution that best fits the distribution of normal- 
ized distances has a = 0.90. This is an indication of a devia- 
tion of parent distributions from Gaussian. 

We can make several conclusions from this test. First, on 
average reported formal uncertainties are correct within sev- 
eral percents. Second, the effect of the core-shift is too small 
to contribute significantly on resul ts of single-epoch surveys. 
According to lKovalev et al.l (l2008h . the typical apparent core- 
shift is expec ted to be 0.4 mas between S- and X-bands. 
iPorcasI (l2009h stressed that when the core-shift is proportional 
to , a source position derived from ionosphere-free linear 
combinations of X and S-band group delays is not sensitive 
to the core-shift and corresponds to a true position of the jet 
base. The core-shift dependence is expected for a conical 
jet with synchrotron self-absorption in the regime of equipar- 
tition between jet particle and magnetic field energy densities 



(lLobanovlll998h . If we assume that this is indeed the case 
for majority of the sources (see also lSokolovskv et al.ll20I lb . 
the average core-shift between K-band and effective S/X posi- 
tions is reduced to 0.4 • , f^/" , , = 0.06 mas. Our observations 

,/k (A -A) 

would allow detecting the core-shift between positions from 
S/X and K-bands observables the 95% confidence level of a 
sample of 190 objects if the variance of the core-shift has been 
greater than 1 mas. VGaPS observations set the upper limit of 
this variance to 1 mas, which does not contradict to results of 
core-shift measurements and predictions. 

9. SUMMARY 

In the VLBA Galactic Plane Survey we detected 327 com- 
pact radio sources not previously observed with VLBI at 
24 GHz in absolute astrometry mode. Half of them are within 
5° of the Galactic plane; 206 of them were also observed and 
detected within the VCS or RDV programs at S/X bands in 
absolute astrometry mode. We determined K-band positions 
of all detected sources. The position uncertainties for all but 
one source are in the range from 0.2 to 20 mas with the median 
value of 0.9 mas. The quoted uncertainties account for various 
systematic effects and their validity within several percents 
was confirmed by comparison with independent X/S observa- 
tions. The detection limit of our observations was in the range 
of 70-80 mJy. For the majority of detected sources parsec- 
scale images were produced and correlated parsec-scale flux 
density estimated. These results are presented in the form of 
the position catalog, calibrated image and visibility data in 
FITS format, and visual plots. 

The new wide-band fringe search baseline-oriented algo- 
rithm for processing correlator output has been developed and 
implemented in the software VIAiA . It has reduces the de- 
tection limit of observations by a factor of \/N, where is the 
number of IF's, by determining group delays, fringe phases at 
the reference frequency, and phase delay rates from the coher- 
ent sum of the data from all IF's. The new algorithm increased 
the number of detected target sources by a factor of 2.4 for this 
survey of weak objects near the Galactic plane. We validated 
the new algorithm by parallel processing of 1080 hours, over 
0.6 million observations, using both the traditional AIPS ap- 
proach and the new approach. The differences between source 
position estimates processed with the wide-field and with the 
traditional AIPS algorithms do not exceed 0.15 mas, which is 
satisfactory for any practical application. 

We investigated possible systematic errors caused by errors 
in the ionosphere path delay derived from GPS TEC maps. 
We derived an empirical model of the ionosphere driven de- 
lay path errors. We found that for declinations > -20° for 
90% of the sources, mismodeling path delay caused source 
position errors of less than 0.05 mas. At declinations below 
-20° these errors grow to 0.15 mas and for some sources may 
reach 0.5 mas. 

Comparison of new K band VLBI positions with positions 
of 192 sources observed at S/X showed an agreement with 
the wrms of 0.46 and 0.6 mas in right ascensions and decli- 
nation respectively, within reported position uncertainties for 
all but two compact steep spectrum sources J0432-I-4138 and 
J2020-I-2942. For these two objects, positional differences 
are about 40 mas. We showed that the reason of these dif- 
ferences is that for sources with complex extended structure 
positions were referred to different source details. These two 
objects demonstrate an existence of an overlooked source of 
errors in VLBI position catalogues that will be studied in de- 
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tail in the future. A 1 mas upper limit on an apparent core- 
shift effect between 8 and 24 GHz is found for the studied 
sample, in agreement with core-s h ift measurements and pre - 
dictions bv Ikovalev et alJ (l2008h . ISokolovsky et alJ (|2011l) . 
iPorcasI (120091k 
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